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Summary 

During  the  three  project  years,  empirical  studies  of  scale  relationships  in  the  retrieval  of  sea  ice 
lead  statistics  have  been  undertaken,  as  have  modeling  investigations  of  atmospheric  influences  on  the 
satellite  signal.  Additionally,  we  have  developed  statistical  models  that  describe  the  scaling  properties  of 
leads.  The  empirical  studies  have  been  based  primarily  on  comparisons  within  and  between  Landsat  and 
AVHRR  imagery,  while  the  atmospheric  models  have  been  specific  to  the  AVHRR.  Submarine  sonar  data 
have  been  used  in  the  statistical  model  development. 

Specific  accomplishments  include:  atmospheric  temperature  and  humidity  profiles  for  the  Arctic 
have  been  constructed  from  Soviet  ice  island  data  and  were  used  in  the  construction  of  three-season 
"standard"  atmospheres  for  the  central  Arctic;  resampling  methods  have  been  tested  on  various  types  of 
imagery,  and  nearest-neighbor  resampling  has  been  shown  to  be  the  most  effective  in  maintaining  the 
spectral  characteristics  of  leads  while  spatial  interpolation  (e.g.,  bilinear)  retains  their  spatial  structure; 
empirical  relationships  between  pixel  size  and  lead  width  have  been  determined;  procedures  for  the 
retrieval  of  lead  statistics  have  been  developed  and  applied  to  Landsat  and  ERS-1  SAR  data;  the 
relationship  between  "apparent"  lead  widths  measured  along  a  transect  (e.g.,  from  submarine  sonar  or  as 
a  sampling  method  for  satellite  imagery)  and  the  "true"  lead  width  distribution  has  been  formalized  in  a 
statistical  sense,  so  that  one  distribution  may  be  obtained  from  the  other;  a  statistical  model  has  been 
developed  for  the  retrieval  of  lead  area  fraction  from  measurements  along  a  line;  a  method  to  partially 
adjust  the  total  lead  area  fraction  for  field-of-view  was  developed  and  tested  on  AVHRR  and  Landsat 
imagery;  a  combination  of  SAR  and  AVHRR  was  used  to  estimate  true  lead  temperature  and  width;  the 
effects  of  atmosphere/surface  conditions  on  the  AVHRR-measured  radiance  in  the  thermal  channels  have 
been  examined  in  terms  of  the  thermal  contrast  between  leads  and  the  surrounding  ice  pack  and  the 
relationship  between  atmospheric  optical  depth  and  lead  size  has  been  quantified;  information  on  Arctic 
aerosol  optical  depths  during  LEADEX  has  been  collected  by  two  of  the  investigators  and  has  been  used 
to  assess  the  extent  of  tropospheric  and  stratospheric  aerosols  and  their  effect  on  lead  detection  from 
satellite;  satellite  and  in  situ  data  collected  during  LEADEX  and  SIMMS ’92  have  been  used  to  test  the 
theoretical  and  empirical  models/methods  developed  during  the  first  two  project  years,  and  shows  these 
models/methods  to  be  generally  valid. 

Additionally,  two  workshops  for  the  satellite  remote  sensing  investigators  of  the  Leads  ARI  were 
hosted  by  this  group  in  Boulder.  One  workshop  report  and  seven  referreed  papers  have  been  published, 
with  three  others  in  press  or  submitted  for  publication.  Two  graduate  students  have  been  supported  part- 
time  over  the  course  of  the  project. 
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Project  objectives 

The  goal  of  this  project  is  to  understand  how  sensor  characteristics,  atmospheric  properties,  and 
surface  conditions  influence  the  detection  and  interpretation  of  sea  ice  leads  using  Advanced  Very  High 
Resolution  Radiometer  (AVHRR)  and  other  satellite  data.  We  seek  to  determine  the  sources  and 
magnitudes  of  errors  inherent  in  the  measurements,  how  data  from  different  sensors  can  be  combined,  and 
how  lead  statistics  change  with  the  different  spatial  resolutions  of  existing  and  future  sensors. 

In  partial  fulfillment  of  these  objectives,  we  have  defined  which  atmospheric  and  surface 
parameters  are  most  critical  for  lead  detection.  Based  on  model  simulations,  we  have  been  able  to  better 
evaluate  the  importance  of  Arctic-specific  model  parameters;  e.g.,  temperature  and  humidity  profiles  and 
aerosols,  for  the  purpose  of  accurately  simulating  sensor  responses.  Sensor  characteristics  such  as  spectral 
response,  field-of-view,  spatial  resolution,  scan  geometry  and  data  processing  methods  coupled  with  scene 
variability  (solar  zenith  angle,  atmospheric  opacity,  surface  temperature,  snow  cover,  ice  thickness,  size 
with  respect  to  sensor  resolution)  determine  feature  signatures.  The  effects  of  these  parameters  had  to  be 
examined  before  lead  signatures  could  be  evaluated  in  terms  of  lead  width  and  orientation,  particularly 
for  features  that  occur  near  spatial  and  radiometric  limits  of  sensor  resolution. 

The  NOAA  AVHRR  satellite  sensor  provides  daily,  Arctic-wide  coverage  of  ice  conditions  at 
moderate  resolution  and  low  cost.  These  image  sets  contain  information  that  is  of  primary  concern  to 
research  and  operational  interests  in  the  Arctic.  Although  a  variety  of  studies  have  examined  various 
aspects  of  remote  sensing  of  sea  ice,  essentially  no  work  had  previously  been  done  to  relate  lead  signatures 
observed  in  AVHRR  data  to  lead  characteristics.  This  lack  of  substantive  verification  work  left  key 
questions  unanswered  and  posed  significant  research  problems  relevant  to  current  lead  investigations. 
Specifically,  the  following  questions  motivated  our  research: 

Lead  and  Surface  Characteristics-.  How  does  lead  detection  depend  on  ice  thickness  for  given  sets 
of  sensor  properties,  surface  temperatures,  and  atmospheric  conditions?  Since  the  temperature 
contrast  between  open  water  and  ice  provides  a  means  to  map  leads  using  thermal  imagery,  to 
what  degree  does  this  contrast  affect  the  apparent  width  of  a  lead  as  observed  in  an  image  and  our 
ability  to  detect  it?  Does  lead  orientation  affect  lead  detection  when  a  wide-angle  scanning 
instrument  such  as  AVHRR  is  used  instead  of  a  nadir-viewing  sensor  such  as  Landsat?  How 
accurately  must  surface  temperatures  be  measured  to  yield  accurate  lead  calculations? 

Atmospheric ,  Boundary  Layer,  and  Solar  Zenith  Angle  Effects :  How  do  these  factors  combine 
with  surface  conditions  and  path  length  to  the  sensor  to  determine  the  thresholds  of  lead  detection? 
What  feedbacks  to  the  atmosphere  do  leads  create,  and  how  will  these  affect  detection;  e.g.,  ice 
crystal  plumes  from  open  leads  that  extend  up  to  -  and  in  some  cases  through  -  the  top  of  the 
inversion  layer?  What  are  the  characteristics  of  "typical"  polar  atmospheres  (i.e.,  water  vapor 
content,  temperature  profiles,  cloud  microphysical  characteristics),  how  are  they  treated  in  radiative 
transfer  models,  and  how  do  these  factors  affect  remote  sensing? 

General  Sensor  Considerations :  In  what  ways  might  sensor  scan-angle,  sensor  calibration,  data 
gridding,  and  image  enhancements  influence  the  ability  to  detect  leads?  Are  lead  statistics  derived 
from  image  centers  where  spatial  resolution  is  greatest  comparable  to  those  derived  from  image 
limbs  where  resolution  is  poorest? 

While  we  do  not  claim  to  be  able  to  answer  all  of  these  questions  completely,  the  results  of  the  past  three 
years’  research  have  given  us  at  ieast  partial  answers  to  each. 
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To  accomplish  our  objectives,  our  approach  included  both  modeling  and  empirical  studies  (Figure 
1.1).  The  radiative  transfer  modeling  is  done  for  the  purpose  of  simulating  the  satellite  sensor  response 
under  a  variety  of  atmospheric  and  surface  conditions.  The  empirical  studies  include:  a  comparison  of 
lead  statistics  determined  in  imagery  of  varying  scales;  e.g.,  AVHRR,  Landsat,  SAR,  and  OLS  data;  effects 
of  different  resampling  methods  and  digital  image  enhancements  on  lead  detection  in  AVHRR  and  Landsat 
data;  use  of  distributions  derived  from  the  low  resolution  imagery  to  estimate  characteristics  of  the 
distributions  obtained  in  the  high  resolution  images;  and  the  relationship  between  lead  width  and  spacing 
statistics  measured  along  a  transect  to  the  true  distributions. 

Results  from  this  woik  will  be  important  to  the  development  and  application  of  lead  detection  and 
mapping  algorithms  proposed  elsewhere  within  the  Leads  ARI.  For  example,  the  ability  to  more 
accurately  access  lead  width  and  spacing  distributions  from  medium  resolution  imagery  is  crucial  to  the 
evaluation  of  large-scale  heat  flux  estimates.  The  modeling  and  empirical  approaches  to  quantifying  the 
relationships  of  scale  discussed  here  are  a  necessary  first  step  to  operational  lead  analysis  from  satellite 
data. 

Satellite  data  used  in  this  study  includes  Landsat  visible,  AVHRR  visible  and  thermal,  ERS-1 
SAR,  and  OLS  visible  and  thermal.  Additionally,  lead  statistics  have  been  derived  from  submarine  sonar 
data.  Radiosonde  temperature  and  humidity  profiles  from  arctic  ice  islands  are  employed  for  radiative 
transfer  studies.  Each  data  type  is  described  in  more  detail  in  the  appropriate  sections. 

This  report  is  divided  into  four  parts.  The  first  details  those  studies  that  directly  relate  to 
retrieving  lead  statistics  from  satellite  imagery:  image  preprocessing,  sensor  field-of-view,  and  sampling 
methods.  The  second  part  describes  radiative  transfer  studies  of  how  the  atmosphere  affects  the  satellite 
signal  and  how  this  in  turn  might  affect  lead  statistics  derived  from  the  imagery.  Part  III  details  the  use 
of  data  collected  during  LEADEX  and  other  field  experiments.  Part  IV  summarizes  the  accomplishments 
to  date. 


Fig.  1.1.  Overview  of  the  modeling  and  empirical  approaches  to  the  study  of  lead  mapping  and  relationships  of  scale 
in  satellite  data.  The  radiative  transfer  model  on  the  left  is  used  for  sensitivity  studies  and  to  generate  simulated 
imagery.  This  imagery  as  well  as  Landsat  data  are  degraded  to  coarser  resolutions.  The  simulated  multi -resolution 
data  sets,  along  with  actual  data  from  different  satellite  sensors,  are  used  to  study  the  effects  of  measurement  scales 
on  derived  lead  statistics. 
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PART  I:  IMAGE-RELATED  STUDIES 


The  empirical  portion  of  this  study  involves  the  comparison  of  AVHRR,  Landsat,  and  -  to  a 
limited  extent  -  KRMS  data.  These  sensors  provide  a  broad  range  of  spectral  and  spatial  resolutions.  In 
this  section,  the  potential  effect  of  image  preprocessing  methods  on  lead  statistics  is  examined. 

1.1  IMAGE  PREPROCESSING  AND  ITS  EFFECT  ON  LEAD  STATISTICS 

Figures  1.2  and  1.3  show  the  effect  on  lead  detection  in  AVHRR  imagery  of  different  resampling 
methods  and  different  digital  number  thresholds  for  "lead  /  no  lead"  mapping.  The  plots  in  Figure  1.2  of 
digital  numbers  along  transects  in  AVHRR  data  processed  using  different  resampling  schemes  show  some 
relatively  small  shifts  in  digital  number  (DN)  within  leads  and  adjacent  to  leads,  which  would  yield  a 
change  in  estimated  lead-covered  area  (such  as  at  transect  location  292-293  in  Figure  I.2a).  For  small 
leads,  the  potential  exists  to  perhaps  mask  the  lead  completely  (as  at  transect  location  279  in  Figure  I.2a). 
Figure  I.3a  shows  results  using  a  threshold  that  detects  small  and  large  leads.  The  threshold  used  for 
Figure  I.3b  detects  only  the  larger  leads.  The  effect  of  the  "smoothing"  interpolations  (bilinear  and  cubic 
convolution)  is  to  eliminate  about  40%  of  the  smallest  leads  using  this  threshold  detection  scheme.  Figure 
1.4  summarizes  the  effect  of  threshold  and  resampling  scheme  on  total  lead-covered  area  (in  this  case,  the 
percent  of  an  AVHRR  image  covered  by  leads).  The  effect  of  using  different  resampling  schemes  is  small 
compared  to  the  effect  of  choosing  different  thresholds.  Thus,  vvhile  Figure  1.3  suggests  a  substantial 
reduction  in  small  leads  when  interpolations  are  used  vs.  nearest  neighbor,  the  effect  of  the  loss  of  these 
small  leads  on  total  lead  area  is  relatively  small.  However,  it  is  also  worth  noting  that,  since  the  turbulent 
fluxes  from  a  given  lead-covered  area  may  increase  as  the  proportion  of  lead-covered  area  in  narrow  leads 
(as  opposed  to  wide  leads)  increases,  any  systematic  shift  by  the  processing  scheme  toward  a  reduction 
in  small  leads  may  need  to  be  considered  when  calculating  regional  estimates  of  surface  fluxes. 

To  further  examine  the  effects  of  different  interpolation  schemes,  "synthetic"  images  containing 
a  simulated  lead  or  lead  complex  with  different  shapes  and  dimensions  (e.g.,  the  patterns  shown  in  Figure 
1.5)  are  used.  Four  lead  types  were  created  (e.g.,  "Type  1",  etc.  in  Figure  1.5):  Type  1  represents  narrow 
open-water  leads;  Type  2  wider  open-water  leads;  Type  3  leads  with  new  and  young  ice  growing  from 
the  edges  of  the  leads;  and  Type  4  representing  leads  with  new  ice  building  up  along  one  side  of  the  leads. 
These  images  were  then  resampled  using  bilinear  interpolation  and  nearest-neighbor  methods.  Resampling 
was  applied  to  a  45°  rotation  of  the  original  images.  Table  1.1  lists  the  change  in  lead-covered  area 
resulting  from  the  resampling  and  determined  using  two  different  thresholds.  A  DN  of  10  (a  typical 
reflectance  for  open  water)  was  chosen  to  represent  open  water  and  a  DN  of  85  was  used  for  ice-covered 
pixels.  The  effects  are  clearly  dependent  on  the  threshold  chosen  to  define  the  cut-off  between  whether 
to  consider  a  pixel  as  part  of  a  lead  or  not.  Both  the  bilinear  and  nearest  neighbor  interpolations  cause 
large  reductions  in  the  area  assumed  to  be  open  water.  When  a  threshold  is  chosen  at  a  higher  reflectance, 
such  as  would  be  appropriate  to  detect  all  pixels  with  some  open  water  or  thin  ice  in  them,  then  the 
percentage  of  pixels  with  some  lead-covered  area  increases  substantially  when  bilinear  interpolation  is 
used. 

In  addition  to  affecting  the  calculation  of  lead-covered  area,  the  different  resampling  schemes 
affect  the  appearance  of  ieads  in  an  image,  and  thus  the  ability  of  an  interpreter  or  automated  pattern- 
recognition  scheme  to  detect  the  leads.  In  particular,  nearest  neighbor  resampling  tends  to  break  up  the 
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linearity  of  leads.  Figure  1.6  shows  a  lead  complex  in  Landsat  imagery  resampled  using  nearest  neighbor 
(Figure  I.6a),  with  an  attempt  at  reconstituting  the  linearity  using  a  median  filter  (Figure  I.6b).  As  with 
most  enhancements,  the  improvement  is  subjective,  but  the  lead  patterns  appear  more  well  defined  in  the 
median-filtered  image. 

Additional  insight  into  the  effects  of  spatial  resolution  and  sensor  properties  can  be  gained  by 
comparing  colocated  imagery  from  different  sensors.  Figure  1.7  shows  registered  Landsat  Multi-spectral 
Scanner  (MSS)  imagery  (top)  and  AVHRR  data  (bottom)  for  a  portion  of  the  Beaufort  Sea.  While  the 
same  general  lead  structure  is  apparent  in  both  images,  the  ability  to  detect  the  smaller  leads  is 
considerably  reduced  in  the  AVHRR  image  (maximum  spatial  resolution  of  1.1  km)  vs.  the  80  m  Landsat 
image.  This  effect  of  spatial  resolution  on  lead  detection,  as  well  as  the  effect  of  different  spectral 
information  on  lead  mapping,  is  also  illustrated  in  Figure  1.8,  which  shows  a  subsection  of  the  Landsat 
and  AVHRR  images,  as  well  as  a  KRMS  strip  superimposed  on  the  Landsat  data.  This  comparison  of 
how  leads  are  represented  in  visible-band  wavelengths  (the  Landsat),  thermal  (AVHRR),  and  passive 
microwave  (KRMS),  points  out  the  problems  of  intercomparing  lead  statistics  derived  from  different 
sensors.  In  this  example,  threshold  detection  of  lead-covered  area  in  the  three  data  types  yields  1.1%  lead- 
covered  area  in  the  Landsat,  12.8%  in  the  AVHRR  (which  includes  apparent  low  cloud  with  substantially 
warmer  temperatures  than  the  ice  surface),  and  5.4%  lead-covered  area  in  the  KRMS  image. 

These  different  representations  are  perhaps  better  represented  by  comparing  transects  through  the 
imagery.  Figure  1.9  shows  such  a  transect.  The  lead  located  at  transect  location  32  is  marked  with  an 
arrow  on  the  Landsat  image  in  Figure  1.8.  The  transect  runs  vertically  through  the  imagery.  In  this 
example,  a  contrast  stretch  was  applied  to  the  AVHRR  data  to  enhance  the  subtle  DN  differences  it.  the 
image.  From  these  comparisons  and  a  similar  comparison  of  transects  in  AVHRR  (unenhanceu)  and 
Landsat  (Figure  1.7),  it  is  fairly  clear  that  the  AVHRR  and  Landsat  reveal  similar  lead  p  uterus,  bi  t  that 
the  number  of  leads  detected,  and  the  image  area  that  is  considered  to  be  partially  lead-covered,  is  quite 
sensitive  to  the  DN  threshold  chosen  to  define  lead  area.  In  Figure  1.9,  for  example,  the  large  lead  at 
location  32  could  be  defined  as  having  a  lead  width  from  800  m  to  about  1500  m  in  the  Landsat  data 
depending  on  the  DN  threshold  used,  and  from  300  m  to  1800  m  in  the  AVHRR,  again  depending  on 
which  DN  is  selected.  The  lead  information  contained  in  the  KRMS  data  clearly  is  quite  different  from 
that  shown  in  the  other  image  types. 
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Table  1.1.  Effect  of  resampling  method  (bilinear  interpolation  =  BL,  nearest  neighbor  =  NN)  and  digital 
number  (DN)  threshold  on  estimation  of  percent  lead-covered  area. 


Fractional  Area  Lead  Coverage  (%) 
(by  Resampling  Method) 


Lead  Type 

Threshold 

None 

Bilinear 

NN 

1 

<;  84 

0.65 

1.23 

0.82 

2 

<;  84 

1.83 

2.61 

1.77 

3 

<;  84 

1.83 

2.59 

1.77 

4 

£  84 

1.83 

2.60 

1.77 

1 

<;  io 

0.65 

0.01 

0.00 

2 

<;  io 

1.83 

.19 

0.00 

3 

5  10 

1.83 

0.00 

0.00 

4 

<;  io 

1.83 

0.64 

0.00 

Fig.  1.2.  Change  in  AVHRR  digital  number  (DN)  using  different  resampling  schemes.  Plotted  data  are  two  transects 
(graphs  A  and  B)  through  three  colocated  images  processed  using  nearest  neighbor  (NN)  (solid  line),  cubic 
convolution  (CQ  (dashed  line),  and  bilinear  interpolation  (BI)  (dot-dash  line). 
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Fig.  1.3.  Comparison  of  number  of  leads  encountered  in  AVHRR  imagery  resampled  using  NN,  CC,  and  BI  with 
two  different  DN  thresholds  (graphs  A  and  B)  used  to  represent  leads. 


Kg.  1.4.  Change  in  percent  lead  area  estimated  in  an  AVHRR  image  as  a  function  of  resampling  scheme  and  DN 
threshold. 


TYPE  1  TYPE  3  TYPE  3  TYPE  4 


Fig.  1.5.  Synthetic  images  used  to  assess  effects  of  resampling  on  different  lead  types. 


Fig.  1.7.  Colocated  Landsat  and  AVHRR  imagery  for  the  Beaufort  Sea,  showing  the  representation  of  the  same  lead 
complex  in  the  different  image  types. 


01GITAL  HUMBER 


Fig.  1.8.  Colocated  Landsat,  AVHRR,  and  KRMS  +  Landsat  imagery  for  the  Beaufort  Sea.  The  black  arrowhead 
in  the  Landsat  image  marks  the  location  of  the  lead  and  transect  presented  in  Figure  1.9. 


Fig.  1.9.  Transect  through  the  imagery  in  Figure  1.8. 
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1.2  THEORETICAL  ANALYSIS  OF  THE  EFFECT  OF  SENSOR  RESOLUTION 
ON  OBSERVED  FRACTIONAL  AREA  COVERAGE 

An  analytical  description  of  the  relationship  between  the  satellite-derived  fractional  area  coverage 
of  a  geophysical  field  and  sensor  resolution  is  needed  in  order  to  assess  the  potential  error  in  many 
satellite-derived  products  and  to  understand  in  a  more  quantitative  manner  the  benefits  of  different  sensor 
systems.  While  there  have  been  studies  of  the  effect  of  sensor  resolution  on  parameter  retrieval,  the 
approaches  have  been  empirical  and  specific  to  a  single  geophysical  variable.  In  the  analysis  of  land  cover 
classes,  for  example,  the  variance  within  satellite  images  has  been  examined  as  a  function  of  measurement 
scale  for  the  purpose  of  determining  the  optimal  resolution  for  change  monitoring  (e.g..  Woodcock  and 
Strahler,  1987;  Townshend  and  Justice,  1988).  In  studies  of  cloud  fraction,  real  and  synthetic  data 
containing  cloud  fields  were  degraded  in  resolution,  and  the  fractional  coverage  was  observed  as  a  function 
of  scale  (e.g.,  Wielicki  and  Parker,  1992;  Wielicki  and  Welch,  1986;  Shenk  and  Salotnonson,  1972). 

Even  though  these  studies  are  useful,  no  concise  statement  of  the  relationship  between  fractional 
coverage  and  sensor  resolution  has  been  given,  so  that  the  results  are  difficult  to  generalize  to  other 
geophysical  fields.  A  complete  analytical  description  of  the  problem  is  difficult  at  best,  involving 
geometrical  (viewing  geometry),  spectral  (band  location  and  width),  radiometric  (signal-to-noise  ratio, 
quantization  levels),  and  spatial  (sensor  resolution  or  pixel  size)  properties.  In  this  section  a  first  attempt 
at  an  analytical  approach  to  the  problem  is  described.  We  are  concerned  only  with  the  fraction  of  the 
image  area  covered  by  some  geophysical  parameter;  e.g.,  open  water  fraction.  We  take  an  approach 
similar  to  that  of  Shenk  and  Solomonson  (1972)  where  cloud  fields  were  simulated  and  the  relationship 
between  pixel  size,  cloud  size,  estimated  area  fraction  and  true  area  fraction  were  expressed  for  different 
thresholding  operations.  That  work  is  extended,  however,  by  generalizing  the  problem  to  any  geophysical 
variable  whose  spatial  structure  can  be  described  by  its  autocovariance  function.  Additionally,  a  specific 
probability  density  function  is  used  as  a  model  of  the  subpixel  area  fraction  so  that  the  results  do  not 
depend  upon  a  given  simulation.  In  this  manner  the  results  are  applicable  to  a  wide  variety  of  geophysical 
fields  including  clouds,  sea  ice  fractures  ("leads"),  and  land  cover  classes.  In  the  next  section  the 
analytical  approach  is  described.  The  models  are  then  applied  to  simulated  fields  of  clouds  and  sea  ice 
leads  in  a  satellite  image  context. 

1.2.1  ANALYTICAL  APPROACH 

Our  goal  is  to  determine  the  proportion  of  pixels  in  an  image  that  have  the  characteristic  of 
interest;  e.g..  the  fraction  that  are  cloudy  or  that  are  sea  ice  leads,  etc.,  given  some  thresholding  operation. 
This  depends  on  the  distribution  of  the  subpixel  area  fraction,  which  is  specified  by  its  shape,  mean  and 
variance.  The  variance  depends  on  the  pixel  size  and  the  spatial  structure  of  the  geophysical  parameter, 
described  by  its  autocovariance  function.  The  formalization  that  follows  can  be  applied  to  virtually  any 
geophysical  parameter  whose  spatial  distribution  can  be  described  in  this  manner. 

Let  <?(x)  be  a  measurable  property  (e.g.,  temperature  or  reflectance)  at  a  point  whose  position  is 
represented  by  the  location  vector  x,  and  define  X,  to  be  any  condition  on  q.  For  example,  X,  might  be  the 
condition  q(x)  <  T,-b,  where  T,  is  the  surface  temperature  and  5  is  some  threshold  amount.  The  indicator 
function  /(x)  in  a  square  region  R  is  equal  to  1  if  q(x)  satisfies  X,  and  0  otherwise.  The  fractional  coverage 
for  which  q  satifies  X,  is  given  by 
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P  = 


A? 


J. « 


dx 


(1) 


where  ,4*  is  the  side  length  of  R  and  is  a  normalizing  factor.  For  the  rest  of  the  discussion  R  is  a  satellite 
image.  The  probability  density  function  (pdf)  of  I  is  // 1)  =  P  and  f,(0)  =  1  -  P. 

Now  let  y)  be  a  measurable  property  of  a  pixel  Z  within  R  centered  at  location  y  (again,  a 
vector).  As  measured  by  the  sensor,  qz  would  be  an  average  over  a  pixel: 

Qz(y)  =  A?  q(x)  dx 

where  Az  is  the  side  length  of  the  pixel.  An  argument  could  be  made  for  using  the  sensor  point  spread 
(transfer)  function,  but  for  simplicity  a  rectangular  spatial  response  is  assumed.  The  fractional  area  of  R 
for  which  qz  satisfies  £  is  an  estimate  of  the  true  fractional  coverage  and  is 

p  -  £  !z( y)  (2) 

R 


where  N  is  the  number  of  pixels  in  R  and  Iz  is  the  indicator  function  for  the  pixel  based  on  qz,  defined 
in  the  same  way  as  is  I  (for  a  point)  based  on  q.  Our  goal  is  to  relate  P’  to  P  over  a  range  of  Az. 

To  determine  P'  analytically  the  probability  density  of  lz  must  be  known.  It  is  not  trivial,  and 
depends  on  P,  £,  pixel  size,  and  the  way  in  which  objects  satisfying  £  are  distributed  in  space.  Since  Iz 
is  a  function  of  qz,  which  in  turn  depends  on  the  fractional  coverage  within  a  pixel  pz,  then  (under  certain 
conditions)  Iz  can  be  expressed  in  terms  of  pz.  For  example,  consider  a  cloud  pattern  where  the  cloud  top 
temperature  Tc  is  everywhere  the  same  and  is  less  than  the  surface  temperture  T,.  Let  £  be  a  thresholding 
operation  such  that 


if  qz  <  T.- 6 
otherwise 


where  6  is  some  threshold  amount.  This  is  equivalent  to 

/  .f1  « P,*P  (3) 

z  [0  otherwise 

The  expression  pz~2.  p  states  that  the  fractional  coverage  within  the  pixel  is  greater  than  some  quantity  p, 
which  in  this  example  has  a  value  such  that 

(1  ~P)TS  +  pTc  <  Ta- 6 


In  reality  there  may  be  a  distribution  of  Tc,  although  we  do  not  address  this  issue  here.  So,  based  on  (3), 
the  probability  density  of  Iz  is 


f,,{  1)  -  Prob (P,2  p) 

1,,( 0)  -  Prob(Pz<p)  =  1-yi) 
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where  Pz  represents  the  random  variable  subpixel  fractional  coverage  (with  specific  realization  pj,  and 
Prob(P2  £  p)  represents  the  probability  that  the  fractional  coverage  within  the  pixel  is  greater  than  some 
quantity  p.  Now,  how  is  Prob(Pz  £  p)  determined? 

For  a  single  pixel  the  fractional  coverage  of  the  geophysical  parameter  is 

pjy)  =  Az  W  ck 


which  can  be  used  as  an  estimate  of  the  true  area  fraction  in  the  image  P.  After  Stoyan  et  al.  (1989, 
pg.  184),  the  expected  value  and  variance  of  Pz  are 


H  =  E(PZ)  -  P 

a2  =  Var (PJ  =  E[PZ~E(PZ)]2 

-  E {[A?  j2  /(x)  dx  -  P[[AZ 2  £  /(x')  dx'  -  P|} 

=  Az  jjz  AcXBx-x'll)  dx  dx' 


(4) 

(5) 


where  k,  is  the  autocovariance  function  for  the  indicator  function.  The  effect  of  pixel  size  on  the 
autocovariance  function  has  been  studied  by  Jupp  et  al.  (1988,  1989),  although  the  autocovariance  function 
in  (5)  does  not  depend  on  pixel  size;  i.e.,  it  refers  to  the  true  underlying  (point)  autocovariance.  The 
expression  (5)  for  the  variance  of  the  subpixel  area  fraction  can  be  expanded  as 

o2  =  V  J^r  mr{A\2n-sm  (6) 

-8Ar\/2  cos[7t/4 +^(/)]  +  2r2cos[2^(r)]}  dr 


where 


r 

0,  0  ZtZA 


-Icos'1 

4 


-r/A) 

t/2-1 


A<r^Aj 2 


as  given  in  Rothrock  and  Thorndike  (1984,  with  a  correction  made  here).  This  applies  to  a  square  pixel 
and  an  isotropic  covariance  function. 

If  a  specific  model  distribution  for  Pz  is  assumed,  with  expected  value  and  variance  as  defined 
above,  then  the  density  of  the  pixel  indicator  function  is  also  known.  Here  we  use  the  Beta  distribution, 
a  two-parameter  density  function  defined  over  the  closed  interval  0  £  p  £  1  that  is  often  used  as  a  model 
for  proportions; 


The  two  parameters  can  be  determined  by  maximum  likelihood  estimation  based  on  the  mean  and  variance 
of  the  subpixel  fractional  coverage  in  (4)  and  (5)  (Falls,  1974); 

The  shape  of  the  distribution  is  related  to  the  size  of  the  pixel  relative  to  the  spatial  structure  (e.g.. 
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Wp)  = 


PT-1(1-P)M 

0, 


r(y^) 

r(Y)r(p)’ 


y,P>0;  0<,p<,1 

elsewhere 


cr 


Y  = 


1  -|l 


wavelength)  of  the  geophysical  parameter.  In  the  limiting  case  with  very  large  pixels  relative  to  the 
wavelength  of  cloud  elements,  for  example,  most  pixels  would  have  a  similar  subpixel  cloud  fraction  and 
the  variance  would  be  very  small.  The  distribution  would  therefore  have  a  single  peak.  On  the  other 
hand,  if  the  pixel  size  is  very  small  then  the  likelihood  of  pixels  being  either  completely  overcast  or 
completely  clear  increases,  the  variance  of  the  subpixel  area  fraction  increases,  and  two  peaks  are 
expected.  This  is  illustrated  in  Figure  1. 10  where  the  beta  distribution  is  shown  for  a  mean  area  fraction 
of  0.2  and  variances  of  0.1,  0.05,  and  0.01. 

The  beta  distribution  has  often  been  used  to  describe  cloud  amount  frequency  distributions  (e.g.. 
Falls,  1974;  Henderson-Sellers  and  McGuffie,  1991;  Jones,  1992).  More  recently  a  similar  distribution, 
the  Burger  distribution,  has  also  been  used  (e.g.,  Henderson-Sellers  and  McGuffie,  1991).  The  Burger 
distribution  is  described  by  the  mean  cloud  amount  and  a  correlation  distance  defined  as  the  separation 
distance  between  pixels  at  which  the  autocorrelation  drops  below  0.99.  Correlation  distance  was  also  used 
by  Gringorten  (1979),  who  developed  models  through  simulations  describing  the  probability  of  a 
meteorological  condition  occupying  some  fraction  of  a  line  or  area. 

Jones  (1992)  presents  a  shape  parameter  that  can  be  used  to  describe  the  beta  distribution.  It  is 
defined  as 


m  i 

S> 0.6  implies  a  U-shaped  distribution,  S»0.6  implies  a  flat  distribution,  and  5<0.6  implies  a  distribution 
with  a  central  peak.  Values  of  the  shape  parameter  are  also  shown  in  Figure  1. 10.  In  the  examples  that 
follow,  and  in  most  satellite  remote  sensing  applications  where  pixel  sizes  are  1  km  or  less,  S> 0.6. 

Now  we  return  to  the  estimate  of  the  total  area  fraction  in  an  image,  P',  which  is  the  proportion 
of  pixels  in  the  image  for  which  the  indicator  function  takes  on  a  value  of  1,  as  defined  in  (2).  Given  the 
distribution  of  subpixel  area  fraction  described  here  by  the  beta  pdf,  P'  is  simply  Prob(Pz  £  p)  or  the 
probability  that  the  subpixel  area  fraction  is  greater  than  some  threshold  value  p.  This  is  the  area  under 
the  curve  to  the  light  of  any  given  value  along  the  horizontal  axis  in  Figure  1. 10. 

Figure  1. 1 1  shows  the  estimated  total  area  fraction  for  four  true  area  fractions  as  a  function  of  the 
subpixel  area  fraction  variance  (along  the  abscissa)  and  the  threshold  value.  The  variance  and  the  true  area 
fraction  together  define  the  distribution  of  subpixel  area  fraction  so  that  a  wide  range  of  spatial  structures 
and  pixel  sizes  is  represented  in  the  three  plots,  independent  of  any  particular  geophysical  field.  For  a 
given  autocovariance  function,  pixel  size  decreases  toward  the  right  in  the  plots.  Note  that  there  is  an 


Probability  Density 
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upper  limit  on  the  variance,  defined  by  the  point  at  which  the  two  parameters  of  the  beta  distribution  are 
equal  to  or  less  than  0.  This  point  is  p-p2-  In  theory,  decreasing  the  pixel  size  fraction.  Results  for  three 
subpixel  area  fraction  thresholds  are  shown. 


Subpixel  Area  Fraction 


Fig.  1. 10.  Three  realizations  of  the  Beta  probability  density  function  for  a  mean  subpixel  fractional  coverage  of  0.2 
and  three  different  variances.  The  shape  parameter  S  is  also  given. 
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Thraahold-O  8 


SubpMVwtane* 


Fig.  1.1 1.  Estimated  total  area  fraction  as  a  function  of  the  subpixel  area  fraction  variance  and  true  area  fraction. 
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1.2.2.  APPLICATION 

Now  we  apply  these  models  to  simulated  satellite  data.  Since  these  models  require  some  a  priori 
knowledge  of  the  field’s  spatial  structure,  they  cannot  be  used  to  assess  the  error  in  the  total  area  fraction 
estimate  from  a  single  image  alone.  Instead,  this  section  illustrates  how  the  error  can  be  assessed  for 
typical  realizations  of  two  very  different  geophysical  variables:  clouds  and  sea  ice  leads. 

A  cloud  field  is  simulated  as  a  distribution  of  disks  whose  center  locations  follow  a  binomial  point 
process  and  whose  diameters  are  approximately  normally  distributed  (in  a  true  Gaussian  distribution, 
negative  diameters  would  be  possible).  One  realization  is  shown  in  Figure  1.12,  where  the  mean  diameter 
is  2000  m.  Sea  ice  leads  are  modeled  as  a  Poisson  line  process.  The  mean  spacing  between  the  lines 
(leads)  is  3000  m  and  their  orientations  are  random.  The  lines  are  assigned  thicknesses  (widths)  following 
the  negative  exponential  density  function: 


W  =  je'*X 

where  w  is  lead  width  and  X  is  the  mean  width.  For  the  simulations  X=2(X)  m.  One  realization  is  shown 
in  Figure  1.13. 

(increasing  the  variance)  beyond  this  point  has  no  effect  on  the  estimate  of  the  total  area  fraction. 

To  examine  the  effects  of  pixel  size  on  the  estimated  fractional  area  coverage,  these  images  were 
degraded  by  simple  averaging  of  2  x  2  pixel  cells.  Four  degradations  were  perfomed.  Initial  pixel  size 
is  50  m;  pixel  size  doubles  with  each  degradation  so  that  the  largest  pixel  examined  is  800  m. 

Exponential  covariance  is  a  reasonable  model  for  many  geophysical  parameters  and  is  used  here: 

k,(i)  «  P(1-P)e^' ,  r,aZ0  (7) 

where  a  describes  the  dependence  of  the  covariance  on  the  separation  distance  r.  Implicit  in  this 
expression  is  that  q{\)  is  isotropic.  Using  (7)  in  (6)  gives 

0s  =  (g) 

-&Arj2  cos(n/4  +!;(*)]  -  2r2cos[2?(f)]}  dr 


Note  that  the  effects  of  area  coverage  and  autocovariance  (e.g.,  the  size  of  objects  and  pixel  size)  separate 
out: 


P(1-F) 


which  is  essentially  Jones’s  (1992)  shape  parameter  S.  This  is  not  strictly  true  for  the  correlation  functions 
of  the  cloud  and  lead  models  employed  here,  but  is  still  useful  for  the  purposes  of  this  paper.  The 
parameter  a  in  (7)  and  (8)  can  be  determined  from  observed  autocovariances  by  rewriting  (7)  in  linear 
form: 
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wm  -  ln[P(1  -PA  -ar  (9) 

and  solving  by  least  squares  regression.  For  the  applications  below,  the  parameter  a  is  determined  for 
three  evenly  spaced,  parallel  transects  in  the  imagery  and  then  averaged. 

Table  1.2  illustrates  the  application  of  the  beta  distribution  and  its  estimated  moments  to  the 
synthetic  data  in  Figures  1.12  and  1.13.  Listed  are  the  pixel  size  relative  to  the  smallest  pixel  (where  the 
smallest  pixel  is  1).  a  determined  from  the  image,  the  "true"  area  fraction  P  determined  from  the  highest 
resolution  image  in  which  each  pixel  is  either  empty  or  completely  full,  the  observed  variance  of  the 
subpixel  fractional  coverage  Pz,  and  the  variance  of  Pr  estimated  by  solving  (8)  numerically.  The 
difference  in  the  a  values  for  the  two  different  fields  reflects  their  spatial  structures  where  the 
(auto)covariance  of  the  synthetic  leads  falls  off  more  rapidly  than  that  of  the  clouds.  The  true  area 
fraction  and  the  observed  and  estimated  variances  in  Table  1.2  for  the  cloud  case  were  used  to  generate 
beta  pdfs  for  comparison  with  the  observed  distribution  of  subpixel  aiea  fraction.  The  results  are  shown 
in  Figure  1.14  as  the  complement  of  the  cummulative  probability  distribution  function.  End  effects  are 
due  to  binning  procedures.  Shown  this  way  it  is  straightforward  to  determine  the  total  fractional  area 
coverage  estimate  P'  of  the  image  for  any  threshold.  For  example,  in  the  top  plot  of  Figure  1.14,  any 
threshold  greater  than  0. 1  (in  theory,  any  value  greater  than  0)  would  yield  the  same  cloud  area  fraction: 
0.25.  With  a  larger  pixel  size,  however,  this  is  not  the  case  (Figure  1.14,  bottom).  A  threshold  of  0.2 
produces  a  total  area  fraction  of  0.33  while  a  threshold  of  0.6  yields  a  total  fraction  of  about  0.2. 


Table  1.2.  Relative  pixel  size,  the  covariance  parameter,  true  total  area  fraction,  and  the  observed  and 
estimated  variances  of  the  subpixel  area  fraction. 


Image 

A 

a 

P 

Observed 

Var(/y 

Estimated 

Var  (/>*)' 

Clouds 

1 

0.05 

0.246 

0.185 

0.178 

16 

0.133 

0.124 

Leads 

1 

0.15 

0.033 

0.032 

0.029 

16 

0.010 

0.011 

1  Estimated  using  (8). 
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Fig.  1.12.  A  simulated  cloud  field,  where  cloud 
centers  follow  a  binomial  point  process  and  cloud 
diameters  are  Gaussian-distributed. 


Fig.  1.13.  A  Poisson  line  process,  used  to  simulate 
linear  openings  in  sea  ice. 
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Fig.  1. 14.  Complement  of  the  cummulative  probability  function  of  subpixel  area  fraction  for  the  cloud  field  in  Figure 
1.12,  using  two  different  relative  pixel  sizes.  Shown  are  the  observed  and  estimated  distributions,  where  the 
estimated  functions  are  based  on  the  true  total  area  fraction  and  the  observed  and  estimated  variances  in  subpixel 
area  fraction  (Table  1.2). 
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1.2.3  DISCUSSION 

Except  in  the  case  of  a  uniform  target  and  background,  as  with  a  single  cloud  deck  over  open 
water,  it  will  generally  not  be  possible  to  determine  the  subpixel  area  fraction  and  hence  to  apply  the 
models  presented  here.  For  example,  suppose  thermal  data  are  being  examined  for  the  puipose  of 
estimating  cloud  fraction.  Even  if  the  background  (land  or  ocean)  temperature  is  known,  a  small  deviation 
from  that  temperature  in  a  pixel  could  be  caused  by  a  small  amount  of  very  cold  cloud  or  a  large  amount 
of  cloud  only  slightly  colder  than  the  background,  or  numerous  other  combinations  of  cloud  amount,  cloud 
optical  thickness,  and  cloud  top  temperature.  The  problem  is  that  a  single  threshold  in  terms  of  subpixel 
area  fraction  Pz  can  translate  into  a  range  of  temperature  thresholds  and  vice  versa.  In  theory  there  is  a 
single  Pz  threshold  that  would  yield  the  same  P‘  over  a  wide  range  of  pixel  sizes.  This  can  be  seen  most 
readily  in  Figures  1. 10  and  1. 14.  At  the  smallest  pixel  size  the  distribution  of  pz  is  bimodal  so  that  a  range 
of  thresholds  would  yield  the  same,  correct  P'\  i.e.,  P'  =  P.  At  a  large  pixel  size  there  may  be  only  one 
correct  threshold,  but  it  is  within  the  range  found  for  the  small  pixel  case.  When  dealing  with  DNs 
(temperatures,  reflectances,  etc.),  however,  this  may  not  be  the  case.  Further  research  is  needed 
concerning  the  effect  of  "regularization",  or  the  averaging  over  the  point  spread  function  of  sensors.  The 
work  of  Jupp  et  al.  (1988,  1989)  is  important  in  this  regard. 

Given  the  pixel  unmixing  problem  when  the  spectral  structure  of  the  field  is  complex,  one  way 
(perhaps  the  only  way  at  present)  to  relate  the  DN  threshold  to  the  subpixel  area  fraction  threshold  is  to 
choose  a  DN  threshold  very  close  to  the  background  value.  This  is  analogous  to  choosing  a  small 
subpixel  area  fraction  threshold,  as  in  the  top  plot  of  Figure  1. 1 1 .  If  the  pixel  size  is  small  enough  relative 
to  the  spatial  structure  of  the  field,  then  P'  will  be  a  good  estimate  of  P.  If  the  pixel  is  not  small,  then 
all  that  can  be  said  is  that  P'  £  P  see  Figure  1.1 1,  top).  How  small  is  small? 

To  more  easily  address  this  and  similar  questions,  Figure  1.15  was  constructed  as  an  aid  in  the 
interpretation  of  (8).  The  figure  shows  the  normalized  variance  as  a  function  of  relative  pixel  size  A  and 
the  slope  of  the  exponential  covariance  (unction  a.  The  normalization  was  done  by  computing  (8)  without 
the  P(l-P)  term  so  as  to  remove  the  dependence  on  P.  Therefore,  to  use  this  figure  the  reported  values 
must  be  multiplied  by  this  term  to  retrieve  the  actual  variance  of  the  subpixel  area  fraction  distribution. 
The  two  important  relationships  in  (8)  that  are  illustrated  in  Figure  1.15  are  that  an  increase  in  the  pixel 
size  or  an  increase  in  the  rate  at  which  the  covariance  drops  off  with  distance  both  result  in  a  decreasing 
variance.  Not  shown  in  Figure  1.15  but  implicit  in  (8)  is  the  relationship  between  P  and  the  variance:  the 
variance  of  Pz  is  maximum  when  P=0.5  for  a  constant  A  and  a. 

Now  if  we  interpret  the  previous  condition  that  a  DN  threshold  close  in  value  to  that  of  the 
background  translates  into  a  subpixel  area  fraction  threshold  of  approximately  0.2  (for  the  purpose  of 
illustration),  then  Figures  1. 11  and  1.15  can  be  used  together  to  answer  "what  if'  questions.  For  example, 
suppose  that  the  geophysical  field  had  an  exponential  autocovariance  with  a=0.4  and  the  true  area  fraction 
is  0.2.  What  would  be  the  estimated  total  area  fraction  P'  with  a  relative  pixel  size  of  2?  The  variance 
based  on  Figure  1.15  would  then  be  approximately  0. 1 1  and  from  Figure  1.1 1  the  estimated  area  fraction 
would  be  about  0.28  with  a  threshold  of  0.2,  0.19  with  a  threshold  of  0.5,  and  0.12  with  a  threshold  of 
0.8.  With  a  relative  pixel  size  of  6  the  variance  is  0.056  and  the  area  fraction  estimates  are  0.36,  0.13, 
and  0.02,  respectively. 

In  theory  a  beta  distribution  that  is  consistent  with  the  observed  autocovariance  function  can  be 
chosen.  The  regression  method  of  estimating  the  rate  of  decay  of  the  autocovariance  described  by  (9)  also 
provides  an  estimate  of  P.  From  these  two  parameters  the  variance  in  (8)  is  then  determined,  thereby 
defining  the  beta  distribution.  However,  the  autocovariance  is  affected  by  pixel  size,  as  described  in  Jupp 
et  al.  (1988,  1989).  Therefore,  the  autocovariance  determined  by  (9)  must  be  translated  into  the  "true" 
(point)  autocovariance  before  the  correct  beta  distribution  can  be  determined. 
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Fig.  1.15.  The  variance  of  the  distribution  of  subpixel  area  fraction  computed  from  (8),  normalized  by  P(l-P),  as 
a  function  of  relative  pixel  size  and  the  slope  of  the  exponential  covariance  function  a.  To  determine  the  actual 
variance,  the  contour  value  must  be  multiplied  by  P(l-P). 


Normalized  Variance 
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1.3  EMPIRICAL  ANALYSIS  OF  THE  EFFECT  OF  SENSOR  RESOLUTION 
ON  LEAD  WIDTH  DISTRIBUTIONS 

The  goal  of  this  section  is  to  begin  to  define  the  sources  and  magnitudes  of  errors  in 
retrieved  lead  statistics  as  a  function  of  the  spatial  resolutions  of  existing  and  future  sensors,  and  to  assess 
the  importance  of  these  errors  in  a  physical  context.  Here  we  compare  lead  statistics  retrieved  from 
satellite  imagery  of  varying  spatial  resolution,  and  we  examine  whether  lead  statistics  derived  from 
"medium"  resolution  imagery  (e.g.,  a  field -of-view  of  500  -  1000  m)  can  be  used  to  estimate 
characteristics  of  lead  distributions  that  would  be  obtained  from  higher-resolution  images  with  a  field-of- 
view  of  about  80  m.  The  effects  of  atmospheric  effects  on  the  retrieval  of  lead  geometries  is  examined 
in  Section  II  and  in  Stone  and  Key  (1993). 

1.3.1  METHODS  AND  DATA 

One  way  to  investigate  the  effects  of  sensor  field-of-view  on  retrieved  lead  statistics  is  to  examine 
data  from  existing  sensors  of  various  spatial  resolutions,  an  example  of  which  is  shown  in  Figure  1.8. 
While  the  same  general  lead  structure  is  apparent  in  the  images,  the  smaller  leads  are  obviously  harder 
to  detect  in  the  AVHRR  image  than  both  the  Landsat  image  and  the  KRMS  data.  In  addition  to  the 
differences  in  spatial  resolution,  each  sensor  is  sampling  different  spectral  characteristics.  For  example, 
relatively  thin  ice  forming  within  leads  exhibits  a  low  albedo  and  relatively  high  physical  temperature,  but 
the  microwave  brightness  temperature  differs  dramatically  from  that  of  open  water.  In  this  example,  the 
problems  inherent  in  comparing  lead  statistics  using  imagery  of  different  spatial  resolutions  and  spectral 
characteristics  is  apparent:  application  of  thresholds  to  the  three  data  types  yields  1.1%  lead-covered  area 
in  the  Landsat  image,  12.8%  in  the  AVHRR  (which  includes  apparent  low  cloud  with  substantially  warmer 
temperatures  than  the  ice  surface),  and  5.4%  lead-covered  area  in  the  KRMS  image.  Depending  on  which 
image  type  is  used,  these  lead-fraction  estimates  would  yield  roughly  an  order  of  magnitude  difference 
in  the  estimate  of  turbulent  heat  transported  into  the  atmosphere  from  the  warmer  ocean. 

While  there  are  advantages  to  comparing  lead  statistics  derived  from  different  types  of  imagery, 
such  a  study  is  complicated  by  different  acquisition  times,  spectral  bands  of  the  various  sensors,  and 
geolocation  problems.  To  alleviate  these  sources  of  uncertainty,  we  choose  to  work  with  images  of  a 
single  data  type  that  are  successively  degraded  in  resolution  by  modeling  the  transfer  function  between 
the  initial  data  and  the  desired  resolution  and  then  sub-sampling.  A  spatial  filter  that  estimates  the  point 
spread  function  of  the  Landsat  sensor  is  applied  following  the  methodology  presented  in  Justice  et  al. 
(1989).  At  each  degradation  cycle,  Gaussian  random  noise  is  added  back  into  the  image  to  reduce  the 
smoothing  effects  of  the  filtering  operation. 

We  start  with  Landsat  MSS  band  4  (0.5-0.6  pm)  scenes  of  the  Beaufort  Sea,  March  1988,  with 
an  initial  FOV  of  80  m  (Figure  1.16).  The  fourth-order  trend  surface  is  removed  from  the  original  grey 
scale  image  ( Eppler  and  Full,  1992)  in  order  to  correct  for  brightness  variations  caused  by  typically  low 
sun  angles  in  the  Arctic.  Images  with  FOVs  of  160,  320,  640,  and  1280  m  are  then  created  using  the 
spatial  filter.  Each  degradation  is  segmented  using  a  threshold  based  on  the  Sobel  operator  edge  detector. 
This  procedure  determines  the  discrete  spatial  gradient  at  each  pixel  in  both  dimensions.  When  its 
histogram  is  compared  with  that  of  the  original  image,  the  point  of  intersection  determines  an  adequate 
cut-off  between  a  lead  and  the  ice.  We  note  that  the  lead/not-lead  decision  is  somewhat  subjective;  linear 
features  from  older,  refrozen  leads,  for  example,  may  or  may  not  be  included  as  leads.  To  differentiate 
between  leads  and  other  low-albedo  features  such  as  shadows  and  isolated  open-water  areas,  valid  lead 
fragments  are  identified  using  tests  based  on  width  and  orientation. 
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Linearity  is  determined  through  correlation/regression  analysis  of  pixels  within  the  candidate 
features.  Lead  widths  are  measured  perpendicular  to  the  regression  line,  at  1  km  intervals,  and  the  slope 
of  the  regression  line  is  the  measure  of  the  feature  orientation. 
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In  addition  to  the  series  of  degraded-resolution  Landsat  images,  "synthetic"  images  are  generated 
as  an  additional  tool  to  study  the  effects  of  spatial  resolution  on  observed  lead  characteristics.  These 
synthetic  images  represent  lead  networks  as  recorded  in  a  thermal  channel;  e.g.,  leads  and  the  surrounding 
ice  are  assigned  different  physical  temperatures.  The  reason  for  using  simulated  leads  is  that  their 
geometrical  characteristics  are  completely  known.  Lead  networks  are  simulated  as  a  Poisson  line  process, 
where  the  mean  spacing  of  lines  (leads)  is  2500  m  and  the  orientations  are  random.  The  lines  are  assigned 
thicknesses  (widths)  following  the  negative  exponential  density  function  with  mean  width  X.  For  the 
simulations,  X  -  '*00  m  based  on  lead  observations  derived  from  submarine  sonar  data  (e.g..  Key  and 
Peckham,  1992).  For  this  stage  of  the  study,  leads  can  consist  of  either  open  water  or  thin  ice  within  the 
surrounding  matrix  of  thick  ice.  Three  ice  thicknesses  are  used:  0  (open  water),  5,  and  15  cm;  the 
surrounding  thick  ice  has  a  thickness  of  2  m.  Corresponding  temperatures  are  271,  256,  248,  and  235  K. 
In  the  simulation,  ice  thicknesses  within  leads  are  assigned  probabilities  consistent  with  ice  thickness 
distributions  reported  by  Maykut  (1982).  One  realization  of  the  Poisson  line  process  is  shown  in  Figure 
1. 17.  The  initial  pixel  size  is  137.5  m,  assigned  so  that  the  pixel  size  after  the  third  degradation  is  1. 1  km, 
the  nominal  FOV  of  the  AVHRR  sensor  at  nadir. 

1.3.2  OBSERVED  CHANGES  IN  LEAD  GEOMETRIES  WITH  FIELD-OF-VIEW 

The  distribution  of  lead  widths  corresponding  to  the  images  in  Figure  1.16  (degraded-resolution 
Landsat  imagery)  is  shown  in  Figure  1.18.  The  disappearance  of  small  leads  due  to  reduction  in  contrast 
and  the  apparent  increase  in  the  relative  frequency  of  large  leads  as  pixel  size  increases  can  readily  be 
seen.  In  this  particular  Landsat  sample,  we  find  that  leads  narrower  than  approximately  250  m  disappear 
as  the  resolution  of  the  Landsat  image  is  degraded  to  320  and  640  m.  However,  the  criteria  for  how  a 
given  lead  will  "grow"  in  width  or  disappear  during  image  degradation  depends  on  contrast  in  reflectance 
of  the  lead  compared  to  that  of  the  surrounding  ice.  For  example,  a  narrow,  open-water  lead  might 
increase  in  apparent  width  while  decreasing  in  contrast  as  pixel  size  increases  during  the  first  degradation. 
However,  in  the  subsequent  degradation,  the  lead  may  "disappear"  as  the  averaging  of  the  sub-resolution 
lead  and  the  surrounding  ice  raises  the  pixel  reflectance  above  a  given  threshold.  A  narrow  refrozen  lead, 
in  comparison,  might  disappear  during  the  first  degradation  since  the  brightness  contrast  between  the  thin 
ice  in  the  lead  (rather  than  open  water  in  the  previous  case)  and  the  surrounding  ice  is  initially  smaller. 

Orientations  of  leads  can  also  be  expected  to  change  if  the  orientations  are  anisotropic  (i.e.,  have 
a  preferred  orientation).  An  illustration  of  this  is  shown  in  Figure  1.19  for  the  Landsat  image  in  Figure 
1.16.  Results  from  other  Landsat  scenes  show  similar  patterns  and  are  therefore  not  shown.  Lead  widths 
and  orientations  from  the  Simula.^!  lead  networks  (e.g.,  Figure  1. 17)  exhibit  similar  dependencies  on  pixel 
size  although  orientations  do  not  change  substantially  as  with  the  real  data  since  the  basic  pattern  is 
isotropic. 

Figure  1.20  shows  the  change  in  mean  lead  width  as  a  function  of  field-of-view  for  six  Landsat 
images.  We  find  that  while  the  manner  in  which  widths  of  individual  leads  changes  is  highly  variable, 
the  mean  lead  width,  averaged  over  the  entire  image,  seems  to  change  in  a  more  predictable  way.  This 
is  a  potentially  important  property  since  "true”  mean  lead  widths  might  then  be  predicted  based  on 
measurements  from  a  lower-resolution  sensor. 

The  change  in  total  lead  areal  coverage  as  a  function  of  field-of-view  is  illustrated  in  Figure  1.21 
for  the  six  Landsat  images.  The  change  in  area  fraction  with  increasing  pixel  size  is  generally  exponential. 
The  actual  rate  of  change  is,  however,  sensitive  to  the  threshold  levels  used.  In  fact,  it  can  be  shown  that 
lead  area  fraction  may  either  increase  or  decrease  with  increasing  pixel  size  depending  on  the  threshold 
used.  However,  when  the  same  thresholding  method  is  used  the  lead  fraction  difference  between 
degradation  cycles  varies  in  a  predictable  way.  In  other  words,  though  the  definition  of  a  "lead”  in  the 
original  image  is  still  subjective,  once  defined  it  remains  consistent  throughout  the  range  of  degradations. 
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Fig.  1.17.  One  realization  of  a  lead  network  simulated  by  a  Poisson  line  process  with  thick  lines.  Pixels  sizes  are 
137.5  m  (upper  left),  275  m  (upper  right),  550  m  (lower  left),  and  1 100  m  (lower  right).  Grey-scale  values  represent 
brightness  temperatures. 
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Fig.  1.18.  Lead  width  distributions  for  the  Landsat  image  series  in  Figure  1.16.  Widths  are  grouped  in  100  m  bins. 
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Fig.  1.19.  Lead  orientations  for  the  degraded  Landsat  series  shown  in  Figure  1.16.  Orientation  is  the  angle  that  a 
lead  makes  with  the  horizontal  axis,  measured  counterclockwise. 
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The  reason  for  the  change  in  lead  geometrical  characteristics  with  sensor  resolution  is  now 
examined  in  terms  of  contrast.  We  define  the  normalized  contrast  as  a  ratio  based  on  the  target  and 
background  temperatures,  Tt  and  TB: 


c-  Irh 

tb 


Of  course,  the  contrast  ratio  need  not  be  defined  in  terms  of  temperature,  so  that  Tr  and  TB  could  also  be 
reflectance  or  digital  number  (ON).  Letting  p  be  the  fractional  area  coverage  of  a  lead  in  a  pixel;  e.g., 
width/FOV,  then  the  total  contrast  which  takes  into  account  the  reduction  in  temperature  contrast  as  a 
function  of  pixel  size  is 

_  [prr+(i-p)re]-re 

^ 'tot  ~  - Y - 
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The  change  in  total  contrast  can  be  seen  in  Figure  1. 16,  and  is  shown  in  more  detail  in  Figure  1.22 
where  four  individual  leads,  each  with  a  different  initial  contrast,  are  placed  in  an  image  context.  The  lead 
at  the  top  of  the  figure  has  the  lowest  initial  contrast.  The  images  are  degraded  as  described  previously, 
with  noise  added  initially  and  at  each  degradation.  The  change  in  the  total  contrast  of  each  lead  from  one 
degradation  to  the  next  is  shown  in  Figure  1.23. 

If  every  pixel  in  the  image  is  to  be  labeled  as  either  a  lead  pixel  or  not  a  lead  pixel,  then  some 
thresholding  operation  must  be  used.  One  possible  method  is  to  choose  as  a  threshold  the  background 
temperature  plus  some  multiple  of  its  variability  a,  say  TB+2a.  This  threshold  can  also  be  expressed  as 
a  unitless  contrast  ratio: 


If  the  total  contrast  of  a  pixel  is  below  this  value,  then  the  pixel  is  not  a  lead  pixel.  This  threshold 
contrast  includes  implicitly  the  effect  of  the  fractional  area  coverage  of  a  lead  within  the  pixel.  It  can  be 
used  to  determine  the  minimum  initial  contrast,  or  critical  contrast,  necessary  for  a  lead  of  a  given  width 
in  a  pixel  of  a  given  size  to  be  detectable: 


P 


where  the  asterisk  represents  a  critical  (cutoff)  value  and  C*lo,=y, 

Figure  1.24  shows  total  contrast  as  a  function  of  the  initial  contrast  and  the  width/FOV  ratio  of 
leads;  i.e.,  the  combinations  of  the  later  two  variables  that  give  rise  to  a  specific  total  contrast.  For 
example,  an  initial  contrast  of  0.15  and  a  p  (width/FOV)  of  0.15  yields  the  same  total  contrast  (0.02)  as 
an  initial  contrast  of  0.05  and  a  p  of  0.4.  The  total  contrast  can  also  be  considered  as  the  threshold 


Part  /:  Image-related  Studies 


35 


contrast  in  that  any  point  below  a  contour  chosen  as  the  threshold  contrast  represents  a  lead  that  is  not 
detectable.  For  example,  given  a  background  ice  temperature  of  240  K  with  a  standard  deviation  of  5  K, 
the  threshold  contrast  as  defined  above  is  0.042.  If  there  exists  a  lead  that  is  500  m  wide  passing  through 
a  1  km  pixel,  (so  p  =  0.5)  then  its  initial  contrast  must  be  at  least  0.084  (which  is  its  critical  contrast)  for 
it  to  be  detected.  Given  a  background  temperature  of  240  K,  this  critical  contrast  translates  into  a  lead 
temperature  of  260.2  K. 


Fig.  1.22.  Four  single-lead  images  of  varying  initial  contrast  (not  shown)  degraded  three  times.  Initial  lead  width 
is  one  pixel.  Gaussian  noise  is  added  after  each  degradation. 


Fig.  1.23.  Total  contrast  as  a  function  of  field-of-view  for  the  four  leads  in  Figure  1.8. 


Total  and  Threshold  Contrast 


Width/FOV 

Fig.  1.24.  Total  contrast  as  a  function  of  the  initial  contrast  and  the  width/FOV  ratio  of  leads;  i.e..  the  combinations 
of  the  later  two  variables  that  give  rise  to  a  specific  total  contrast.  The  total  contrast  can  also  be  treated  as  the 
threshold  contrast  in  that  any  point  below  a  contour  chosen  as  the  threshold  contrast  represents  a  le.  that  is  not 
detectable. 
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1.3.4  EFFECTS  OF  CHANGES  IN  FIELD-OF-VIEW  ON  TURBULENT  HEAT  FLUX 

As  noted  in  previous  sections,  image  characteristics  affect  the  retrieval  of  lead  width  distributions 
as  well  as  total  lead  area.  To  quantify  the  effects  of  these  image  characteristics,  we  calculate  turbulent 
heat  flux  as  a  function  of  lead  area  and  lead  width.  Changes  in  both  the  mean  lead  width  and  lead- 
covered  area  are  considered  in  the  calculation  of  sensible  and  latent  heat  flux  as  a  function  of  fetch 
(treated  here  as  the  lead  width),  surface  temperature,  air  temperature,  and  wind  speed  using  the  procedure 
outlined  by  Andreas  and  Murphy  (1986).  In  this  approach,  a  bulk  Richardson  number  defines  atmospheric 
stability  that  controls  convective  turbulence  based  on  temperature  and  wind  speed.  Convective  turbulence 
combines  with  the  mechanical  mixing  introduced  by  the  step  effect  of  an  air  mass  in  equilibrium  with 
thick  sea-ice  conditions  travelling  over  the  physically  rough  edge  of  a  lead  and  the  considerably  warmer 
open  water  or  thin  ice  in  the  lead.  The  addition  of  mechanical  turbulence  introduced  by  the  ice-lead 
boundary  tends  to  result  in  a  higher  rate  of  heat  transfer  from  smaller  leads  compared  to  larger  leads. 
Thus,  for  a  given  areal  coverage  of  leads  in  an  image,  a  greater  number  of  smaller  leads  will  result  in 
more  heat  loss  to  the  atmosphere  than  from  a  lesser  number  of  larger  leads,  even  though  the  total  amount 
of  open  water  in  the  image  remains  the  same.  Under  the  conditions  examined  by  Andreas  and  Murphy 
( 1986),  this  decrease  in  flux  as  lead  width  increases  becomes  negligible  for  lead  widths  greater  than  about 
200  m. 

To  illustrate  the  effects  of  changes  in  lead  statistics  using  different  image  fields-of-view,  we 
calculate  sensible  and  latent  heat  flux  using  the  above  approach  for  the  data  presented  in  Figure  1.20  and 
the  associated  changes  in  lead  areal  coverage  in  Figure  1.21.  An  open-water  temperature  of -1.8°  C,  wind 
speed  of  5  m  s'1,  air  temperatures  of  -28.9°  C  at  a  reference  height  of  2  m,  ocean  salinity  of  34  ppt,  air 
pressure  of  1000  mb,  and  a  neutral-stability  drag  coefficient  of  1.49  x  10 3  are  used  to  represent  typical 
mid-winter  (January)  conditions  over  the  Arctic  sea  ice  pack  (Maykut,  1978;  Andreas  and  Murphy,  1986). 
Although  leads  are  often  covered  by  thin  ice  rather  than  open  water  and  thus  have  a  lower  surface 
temperature  than  open  water,  the  assumption  that  the  leads  are  not  refrozen  and  have  a  surface  temperature 
of  -1.8°  C  is  a  useful  baseline  for  our  calculations.  Turbulent  (sensible  plus  latent)  heat  flux  from  leads 
is  calculated  using  the  mean  lead  width  at  each  field  of  view  and  then  weighted  by  the  areal  coverage  ~f 
leads  for  the  six  MSS  images  used  in  Figures  1.20  and  1.21  to  yield  an  areally-averaged  heat  flux. 
Turbulent  fluxes  from  open-water  leads  under  these  conditions  are  around  300  Wm 2  compared  to  a  flux 
of  nearly  0  from  surrounding  ice  taken  to  be  three  meters  thick.  Thus,  lead  fraction  and  lead  width 
dominate  the  transfer  of  turbulent  heat  through  the  ice  pack  during  winter.  Table  1.3  shows  these  areal 
averages  for  the  six  MSS  images.  Since  the  effect  of  increasing  the  fields-of-view  in  these  examples  is 
to  decrease  the  apparent  lead  fraction,  areally-averaged  fluxes  decrease  as  field-of-view  increases. 
However,  as  noted  earlier,  the  choice  of  thresholds  can  affects  both  the  magnitude  and  direction  of  change 
in  lead  statistics  with  changing  FOV.  If  we  assume  that  the  lead  widths  and  lead  fractions  measured  using 
the  80  m  FOV  imagery  are  closest  to  reality,  then  the  errors  introduced  by  using  lead  widths  and  lead 
fraction  measured  at  a  640  m  FOV  are  substantial  -  averaging  45%  over  the  six  images.  Since  the  change 
in  turbulent  heat  transfer  with  changing  lead  width  is  greatest  for  smaller  leads,  those  images  with  smallest 
mean  lead  widths  at  the  80  m  FOV  (such  as  Images  C  and  F)  are  most  affected.  In  the  images  studied 
here,  where  the  mean  lead  width  is  fairly  large,  the  effect  of  errors  in  lead  fraction  is  about  five  times  that 
of  the  effect  of  uncertainty  in  lead  width. 
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I 
I 

Table  1.3.  Areally  averaged  turbulent  (sensible  +  latent)  flux  (in  Wm'J)  for  typical  January  conditions  as  m 

a  function  of  field  of  view  for  six  MSS  images.  The  percent  change  in  the  flux  between  FOVs  of  80  m  I 

and  640  m  is  also  shown. 


FOV  (m) 

A 

B 

Image 

C  D 

E 

F 

80 

23.5 

19.1 

17.4 

50.7 

10.6 

15.0 

160 

22.2 

17.7 

14.4 

47.0 

9.7 

12.1 

320 

19.9 

16.1 

11.7 

41.6 

8.2 

8.9 

640 

16.7 

7.9 

33.7 

5.6 

5.3 

%  change 

29 

27 

55 

34 

47 

65 

1.3.5  TRANSLATION  BETWEEN  SCALES 

From  the  previous  discussions  it  is  obvious  that  lead  statistics  change  significantly  as  a  function 
of  field-of-view,  and  that  there  are  important  implications  of  these  changes  for  large-area  turbulent  heat 
flux  estimates.  Is  there  any  possibility  of  estimating  the  true  lead  widths  and  area  fractions  from  those 
observed  in  lower-resolution  imagery? 

1.3.S.1  Width  Distributions 

Given  that  very  small  features  will  generally  not  be  resolved,  the  issue  then  becomes  one 
concerning  the  possibility  of  using  the  distribution  of  lead  widths  measured  at  low  resolution  to  estimate 
the  complete  or  "true"  distribution.  For  example,  assume  that  lead  widths  x  follow  a  negative  exponential 
distribution  with  an  unknown  mean  X.  From  a  sampling  point  of  view  it  is  useful  to  treat  the  distribution 
of  widths  as  discrete  and  address  the  number  n,  of  leads  in  bin  i  that  have  widths  between  x,  and  x,+w: 


n,  = 


NWq  -X,  IX 


(10) 


where  w  is  the  width  of  the  bin  and  N  is  the  unknown  total  number  of  leads  in  the  spatial  area.  The  idea 
is  that  n,  is  measured  for  a  few  bins,  and  that  X  and  N  are  estimated.  To  accomplish  this,  (10)  is  rewritten 
in  linear  form  as 

ln(n)  -  In(^)  -  1*,  (11) 

Letting  a  =  In (NwIX)  and  b  =  A.1  and  solving  for  a  and  b  by  the  method  of  least  squares  with  the  observed 
data,  the  mean  of  the  distribution  and  the  total  number  of  leads  can  then  be  estimated. 

Experiments  with  this  model  show  it  to  be  very  sensitive  to  the  bin  width  and  the  number  of  bins 
in  which  leads  actually  occur  in  the  lower-resolution  imagery.  This  is  not  unexpected  considering  that 
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the  entire  range  of  x  is  being  estimated  in  the  least  squares  model  by  observations  in  only  one  part  of  its 
range.  A  more  fundamental  problem  exists  with  this  method:  as  shown  earlier,  the  widths  of  the  leads 
observed  in  the  lower-resolution  data  are  probably  not  the  true  widths  of  those  leads.  Figure  1.25 
illustrates  the  problem  where  the  actual  lead  width  distribution  -  which  is  exponential  -  in  the  simulated 
leads  image  of  Figure  1.17  is  estimated  using  the  above  method  and  the  obseivations  from  each  degraded 
image.  Significant  departures  from  the  actual  distribution  are  obvious. 


Estimated  Lead  Width  Distributions 


Fig.  1.25.  Estimated  actual  width  distributions  from  observed  lead  widths  at  each  field-of-view  or  degradation 
(DEGO-DEG3)  of  the  simulated  lead  network  shown  in  Figure  1.17.  Also  shown  is  the  model  distribution  used  to 
generate  the  leads  in  the  image. 
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Is  it  possible  to  unmix  the  pixels  and  thereby  obtain  the  true  lead  widths?  Using  a  single  spectral 
band,  it  is  not  possible  when  ice  of  different  thicknesses,  and  thus  different  reflectances  and  surface 
temperatures,  are  present  in  the  field-of-view.  For  a  brief  time  during  the  summer  when  new  ice  is  not 
forming  in  the  leads,  the  percentage  of  open  water  w»thin  the  FOV  can  be  calculated  with  a  single  spectral 
band  since  all  leads  can  be  assumed  to  contain  only  open  water  and  therefore  essentially  the  same 
reflectance  or  temperature.  During  the  winter  when  no  visible-band  data  are  available,  no  unmixing  is 
possible  since  leads  can  consist  of  a  large  range  of  ice  thicknesses.  During  the  spring  and  fall  months  the 
problem  can,  in  theory,  be  solved  using  one  thermal  and  one  visible-band  observation  and  an  energy 
balance  approach  as  follows.  The  total  contrast  of  a  lead  pixel  in  both  a  thermal  and  a  visible-band  image 
of  the  same  lead  are  observed.  The  mean  background  (ice)  temperature  and  albedo,  T0  and  aB  are 
determined  from  the  data.  This  leaves  two  equations  with  three  unknowns: 


^ tot  JR 

0,0t,  Vj3 


=  P 


=  P 


Tt-Tb 

tb 

ar-as 

as 


Actually,  the  target  (lead)  temperature  and  albedo  are  physically  related,  although  the  relationship  is  a 
complex  one.  An  energy  balance  model  is  used  to  determine  the  target  albedo  for  a  given  target 
temperature  (Maykut,  1982): 


(1-cc T)Fr-liCe+FL+e  oT$+Fa+Fg  +  Fc  =  0 

where  a  is  the  albedo,  e  is  the  longwave  emissivity,  o  is  the  Stefan-Boltzmann  constant,  /ice  is  the  amount 
of  shortwave  energy  that  penetrates  the  ice  and  does  not  directly  heat  the  surface,  Fr  and  F,  are  the 
downwelling  shortwave  and  longwave  radiation,  F,  and  Fe  are  the  sensible  and  latent  heat  fluxes,  and  Fr 
is  the  conductive  heat  flux.  A  flux  toward  the  surface  is  positive.  The  energy  balance  equation  is  solved 
for  a  range  of  possible  target  temperatures,  TB  £  Tr  £  273.15  K,  until  a  combination  of  p,  Tt,  and  a,  is 
found  that  is  consistent  with  the  observed  total  contrasts. 

While  in  theory  this  method  will  work,  in  practice  it  would  be  difficult  to  accurately  estimate  all 
the  necessary  parameters.  It  is  not  our  purpose  here  to  present  methods  of  retrieving  these  parameters, 
but  instead  we  summarize  the  potential  error  through  an  example  of  the  sensitivity  of  the  energy  balance 
approach:  if  the  target  albedo  can  be  estimated  to  within  0.05,  for  example,  the  range  of  p  that  could 
satisfy  the  above  equations  is  0.455  to  0.556  for  a^O.7,  a  true  p  of  0.5  and  a  true  of  0.2.  With  a  1  km 
FOV  this  translates  into  a  range  in  lead  widths  of  445  to  556  m,  where  the  true  width  is  500  m.  While 
the  use  of  physical  models  can  help  the  unmixing  process,  we  do  not  expect  that  we  will  ever  be  able  to 
fully  resolve  the  mixture  components  with  existing  data. 

As  can  be  seen  in  Figure  1.20,  however,  the  unmixing  of  pixels  to  determine  the  actual  lead  widths 
observed  may  not  be  necessary.  The  fact  that  the  mean  lead  widths  change  in  a  predictable  way  with 
increasing  pixel  size  implies  that  the  mean  of  the  width  distribution  measured  at  one  field-of-view  can  be 
used  to  estimate  the  mean  width  at  another  field-of-view.  Even  though  the  rate  of  change  of  mean  lead 
width  with  pixel  size  depends  on  the  threshold  used  (not  shown),  the  relationship  is  approximately  linear 
for  a  given  thresholding  operation.  This  relationship  can  be  utilized  as  follows.  For  an  image  at  a  given 
field-of-view,  determine  the  mean  lead  width.  Degrade  the  image  and  once  again  determine  the  mean  lead 
width  using  the  same  thresholding  operation  as  before.  The  two  points  define  a  line  analogous  to  those 
in  Figure  1.20.  The  mean  lead  width  at  a  narrower  FOV  can  then  be  determined  and  applied  in  (11). 

Of  course  the  relationship  is  not  perfectly  linear,  so  that  some  error  in  the  predicted  mean  lead 
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width  can  be  expected.  As  an  example,  we  consider  the  two  images  with  the  smallest  and  greatest  error 
in  turbulent  heat  flux  due  to  FOV,  as  given  in  Table  1.3,  images  B  and  F.  The  relationships  for  these  two 
images  are  among  the  most  nonlinear  of  those  examined.  For  image  B  the  mean  lead  width  at  a  pixel  size 
of  320  m  is  680  m;  degrading  that  image  to  a  640  m  pixel  size  would  yield  a  mean  lead  width  of 
approximately  1080  m,  also  as  indicated.  If  these  two  points  were  then  used  to  extrapolate  back  to  a  FOV 
of  80  m,  the  estimated  lead  width  would  be  360  m  as  opposed  to  the  23 1  m  value  actually  measured  in 
the  imagery.  The  difference  between  the  turbulent  flux  calculated  for  a  single  lead  (no  area-weighting) 
whose  width  is  the  mean  width  at  320  m  FOV  and  that  calculated  at  the  80  m  FOV  is  16  Wm'2  (4.1% 
difference)  compared  to  12  Wm'2  (3.1%)  using  the  mean  width  extrapolated  to  an  80  m  FOV  from 
measurements  at  320  m  FOV.  For  image  F,  which  showed  the  greatest  sensitivity  of  turbulent  fluxes  to 
field-of-view,  the  corresponding  errors  are  25  Wm'2  (6.3%)  and  15  Wm'2  (3.8%). 

13.5.2  Total  Area 

If  lead  area  fraction  in  satellite  imagery  follows  a  known  scaling  law,  then  the  "true"  area  fraction 
can  be  estimated  from  the  area  fraction  determined  at  any  scale.  Fractal  geometry  is  type  of  scaling 
relationship  that  has  been  used  in  the  analysis  of  geophysical  phenomena  and  deserves  mention  in  the 
context  of  lead  area.  In  particular,  the  stream  length-drainage  area  relationship  has  been  described  in  terms 
of  fractals  (Robert  and  Roy,  1990).  However,  in  that  and  related  studies  the  streams  have  no  width  and 
are  therefore  not  applicable  to  the  lead  studies  presented  here.  In  contrast,  Karlinger  and  Troutman  (1992) 
have  examined  the  "fat"  fractal  relationship  between  river  channels  with  finite  widths  and  drainage  area. 
An  examination  of  the  data  presented  in  Figure  1.21  reveals  that  in  general  the  fractional  area  coverage 
of  leads  decreases  exponentially  (log-linear)  with  increasing  pixel  size,  and  in  some  cases  the  decrease  is 
even  linear,  so  that  the  log-log  relationship  described  by  fractal  scaling  laws  does  not  appear  to  apply. 

As  with  lead  widths,  the  rate  of  change  of  area  fraction  with  increasing  pixel  size  is  not  constant, 
but  rather  is  a  function  of  the  threshold  used.  In  fact,  the  direction  of  change  is  also  threshold-dependent, 
so  that  the  lead  area  may  increase  or  decrease  with  increasing  pixel  size.  The  theoretical  reasons  for  this 
are  examined  in  Key  (1993)  where  the  distributions  of  the  subpixel  area  fraction  of  various  geophysical 
fields  with  known  covariance  structures  are  modeled  by  a  Beta  probability  distribution,  and  the  estimated 
total  area  fraction  in  an  image  is  determined  as  a  function  of  threshold.  Unfortunately,  the  relationship 
between  digital  number  (or  temperature  or  reflectance)  and  the  subpixel  area  fraction  can  be  complex,  so 
that  expressing  the  subpixel  area  fraction  threshold  as  a  DN  threshold  is  often  not  possible. 

Therefore,  perhaps  the  best  estimate  of  the  true  area  fraction  of  leads  in  an  image  is  obtained  using 
the  procedure  outlined  earlier  for  mean  lead  width:  degrading  the  image  once,  assuming  an  exponential 
or  possibly  linear  relationship,  and  extrapolating  back  to  a  smaller  FOV.  Of  course,  the  same  thresholding 
operation  must  be  used  for  both  images  (the  Sobel  operator  here).  As  was  done  in  the  previous  section, 
the  potential  error  in  this  method  for  the  Landsat  images  can  be  examined  using  the  data  in  Figure  1.21. 
Using  image  B,  lead  fraction  extrapolated  to  an  80  m  FOV  from  observations  at  an  FOV  of  320  m  is 
0.031  versus  the  observed  fraction  at  80  m  of  0.0378.  Combining  this  error  with  the  error  in  open-water 
turbulent  flux  associated  with  mean  lead  width  as  calculated  in  the  previous  section,  the  error  in  using  the 
extrapolated  lead  fraction  and  lead  width  versus  the  statistics  observed  at  an  80  m  FOV  is  5.7%  of  the 
areally-averaged  turbulent  heat  flux.  If  the  lead  statistics  are  not  adjusted  for  field-of-view,  i.e.,  if  the  lead 
statistics  observed  at  an  FOV  of  320  m  are  used,  then  the  error  increases  to  15.6%.  For  image  F,  the  error 
is  20%  using  the  extrapolated  statistics  compared  to  40.7%  percent  using  the  statistics  at  320  m.  In  these 
two  cases,  extrapolating  the  lead  statistics  reduces  the  average  error  in  turbulent  flux  by  about  57%.  This 
technique  is  further  applied  to  AVHRR  and  Landsat  TM  data  in  Section  III. 2.2. 
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1.4  LINEAL  METHODS  OF  ESTIMATING  LEAD  PARAMETERS: 

WIDTH  DISTRIBUTIONS  MEASURED  ALONG  A  TRANSECT 

Perhaps  the  largest  source  of  high-resolution  data  that  is  potentially  useful  for  lead  statistics  is  that 
collected  by  submarine  sonar  over  the  last  three  decades.  For  ice  draft  information  these  data  are 
invaluable.  But  can  they  also  be  used  for  statistics  of  lead  geometries:  i.e.,  lead  widths  and  spacings? 
Lead  width  and  spacing  statistics  have  been  examined  in  two  sonar  data  sets.  Both  are  in  the  Canada 
Basin  (Figure  1.26);  one  in  August  of  1970  and  one  in  October  of  1978.  Tables  1.4  and  1.5  show  the 
results.  How  realistic  are  these  data?  If  leads  are  conceptualized  as  linear  features  of  some  width,  then 
crossing  the  lead  at  any  angle  other  than  perpendicular  to  its  local  orientation  will  result  in  an  overestimate 
of  the  actual  width. 


Figure  1.26.  Sonar  transect  locations  through  the  central  Canada  Basin  (August  1970)  and  northeastern 
Canada  Basin  (October  1978). 
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Table  1.4.  Lead  widths  (m)  and  standard  deviations  (below)  in  the  QUEENFISH  data  by  region,  maximum 
draft  (cm),  minimum  width  (m).  Also  given  are  the  maximum  lead  widths  encountered.  No  statistics  are 
given  if  fewer  than  20  leads  were  found  in  a  region. 


MAX  DRAFT/ 
MIN  WIDTH 

A 

B 

C 

REGION 

D 

E 

F 

G 

H 

30/3 

26.5 

33.1 

21.3 

26.5 

47.7 

18.5 

23.8 

2.1 

2.9 

2.4 

2.6 

7.1 

2.6 

3.7 

20 

60.5 

66.1 

54.8 

53.5 

68.3 

3.8 

5.8 

5.6 

5.2 

9.6 

50 

88.2 

109.7 

100.4 

5.3 

10.8 

10.4 

200 


Maximum: 

228 

502 

130 

283 

327 

70 

90 

70/3 

32.6 

54.8 

27.5 

32.1 

49.4 

36.4 

34.7 

40.0 

3.3 

3.8 

2.0 

2.8 

7.0 

8.3 

3.8 

3.7 

20 

63.7 

78.1 

55.6 

54.2 

73.7 

55.9 

59.0 

60.4 

6.8 

5.2 

3.9 

4.2 

9.9 

13.9 

6.4 

5.4 

50 

106.1 

129.2 

106.5 

104.2 

129.5 

98.1 

105.5 

200 

Maximum: 

14.3 

8.5 

7.9 

9.5 

19.0 

10.3 

9.1 

885 

526 

255 

294 

374 

537 

257 

227 

100  /  3 

28.9 

40.8 

32.9 

30.9 

50.5 

32.4 

36.4 

40.0 

2.4 

2.5 

2.5 

2.2 

6.5 

6.6 

3.7 

3.6 

20 

59.3 

70.4 

70.2 

53.6 

74.1 

57.9 

62.8 

5.4 

4.1 

5.5 

3.8 

9.1 

5.9 

5.6 

50 

102.4 

126.6 

124.5 

98.6 

123.0 

98.6 

109.6 

2C0 

Maximum: 

11.9 

7.8 

11.2 

8.1 

15.9 

10.0 

9.8 

886 

533 

763 

299 

374 

539 

263 

269 
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Table  1.4,  cont.  Lead  widths  (m)  and  standard  deviations  (below)  in  the 
northeastern  Canada  Basin  data  by  region,  maximum  draft  (cm),  minimum  width 
(m).  Also  given  are  the  maximum  lead  widths  encountered.  No  statistics  are 
given  if  fewer  than  20  leads  were  found  in  a  region. 


MAX 

MIN 

DRAFT/ 

WIDTH 

I 

J 

REGION 

K 

L 

M 

30 

/  3 

116.9 

100.8 

9.6 

4.9 

in  o  _ 

5.4 

1.7 

1U  .4  - 

2.4 

20 

50 


^UU 

Maximum: 

2181 

153 

83 

42 

70/3 

30.8 

14.4 

17.6 

8.7 

15.9 

16.1 

2.0 

3.3 

1.3 

4.2 

20 

132.8 

62.5 

65.8 

50 

200 

Maximum: 

81.4 

9.3 

13.7 

2274 

234 

342 

72 

510 

100  /  3 

28.7 

17.9 

18.8 

9.4 

16.5 

12.1 

2.1 

2.7 

1.1 

4.5 

20 

104.6 

61.6 

68.8 

84.7 

50 

52.9 

7.4 

10.8 

29.9 

200 

Maximum: 

2280 

239 

348 

81 

792 

I 

I 
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Table  1.5.  Lead  spacings  in  the  QUEENFISH  data,  by  region,  maximum  draft  (cm),  minimum  width  (m). 
S=mean  spacing  (m)  with  standard  deviation  below,  N=number  of  leads  per  kilometer.  Also  given  are 
the  maximum  spacings  encountered.  No  statistics  are  given  if  fewer  than  20  leads  were  found  in  a  region. 


MAX  DRAFT/  REGION 

MIN  WIDTH  ABCDEF  GH 


30/3 

S: 

462 . 9 

428.5 

929.9 

39.9 

41.7 

118.9 

N: 

2.01 

2.16 

1.06 

20 

S: 

1240.1 

1014.0 

2821.6 

150.1 

125.6 

601.9 

N: 

0.73 

0.92 

0.31 

50 

S: 

2552.6 

2351.2 

318.7 

503.1 

200 

Maximum: 

N: 

0.36 

0.41 

4893 

5842 

9321 

70/3 

S: 

339.7 

369.8 

364.1 

24.9 

30.2 

34.6 

N: 

2.67 

2.35 

2.56 

20 

S: 

783.2 

568.2 

926.6 

70.6 

50.3 

97.9 

N: 

1.16 

1.53 

1.02 

50 

S: 

1818.8 

1238.0 

2935.7 

190.6 

150.4 

621.9 

200 

Maximum: 

N: 

0.50 

0.73 

0.33 

3807 

4792 

4777 

100  /  3 

S: 

238.1 

213.4 

225.5 

14.8 

13.1 

17.9 

N: 

3.73 

3.93 

3.86 

20 

S: 

597.7 

416.5 

581.6 

46.4 

33.7 

59.9 

N: 

1.51 

2.05 

1.54 

50 

S: 

1468.9 

1064.2 

1449.1 

157.0 

126.1 

261.8 

200 

Maximum: 

N: 

0.61 

0.83 

0.63 

2897 

2808 

3628 

739.4 
82.2 

1.31 

1897.2 

307.5 
0.51 

685 . 9 
155.1 
1.38 

1061.8 

240.6 

0.90 

2190.9 

402.5 

0.45 

2551.0 

759.2 

0.33 

5493.8 

1436.6 

0.17 

8594 

7062 

13523 

24321 

513.0 

41.2 

1.83 

504.7 

102.8 
1.82 

553.1 

108.2 
1.70 

859.6 

100.1 

1.10 

737.7 

61.9 

1.14 

1050.0 

102.5 

0.91 

816.9 

173.1 

1.14 

983.4 

221.0 

0.98 

1775.5 

264.0 

0.54 

1230.3 

131.4 

0.67 

3226.0 

577.5 

0.29 

1868.1 

515.6 

0.48 

3757.9 

1200.8 

0.22 

3029.1 

653.3 

0.27 

4488 

6782 

4996 

6221 

3304 

416.3 

31.5 

2.24 

443.2 

88.4 

2.04 

421 . 6 
84.4 
2.20 

837.0 

97.6 

1.13 

722.0 

72.3 

1.31 

866.6 

78.7 

1.09 

705.2 

147.0 

1.30 

1587.7 

227.5 

0.60 

1321.9 

213.0 

0 . 72 

2465.7 

348.7 

0.37 

1591.7 

397.6 

0.60 

3603.1 

1158.5 

0.23 

3251.1 

723.5 

0.29 

4478 

6779 

4993 

6218 

8918 
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Table  1.5,  cont.  Lead  spacings  in  the  northeastern  Canada  Basin  data,  by 
region,  maximum  draft  (cm),  minimum  width  (m).  S=mean  spacing  (m) 
with  standard  deviation  below,  N=number  of  leads  per  kilometer.  Also 
given  are  the  maximum  spacings  encountered.  No  statistics  are  given  if 
fewer  than  20  leads  were  found  in  a  region. 

MAX  DRAFT/  REGION 


MIN 

WIDTH 

I 

J 

K 

L 

M 

30 

/  3 

S: 

3381.9 

4567.3 

2477.4 

-  2996.1 

1500.2 

2495.0 

913.9 

1948.1 

N: 

0.15 

0.22 

0.36 

0.17 

20 


DU 

Maximum: 

23429 

71342 

42161 

37059 

70/3 

S: 

994.4 

700.1 

928.5 

1763.9 

911.1 

229.6 

129.7 

159.7 

449.3 

169.6 

N: 

0.96 

1.37 

1.01 

0 .53 

1.07 

20 

4476.4 

4587.8 

4479.7 

1396.3 

1513.8 

1513.3 

50 

200 

Maximum: 

N: 

0.19 

0.21 

0.20 

24918 

13508 

10022 

25355 

14562 

100  /  3 

S: 

739.6 

585.6 

635.4 

858.2 

709.5 

130.6 

73.9 

89.7 

197 .8 

110.9 

N: 

1.28 

1.65 

1.53 

1.09 

1.36 

20 

S: 

2896.8 

2765.0 

2860.5 

4493.1 

795.8 

623.6 

608.4 

1214.4 

50 

200 

Maximum: 

N: 

0.29 

0.35 

0.31 

0.19 

12353 

7298 

10022 

25212 

11154 
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A  methodology  has  been  developed  for  estimating  width  distributions  of  linear  features  from 
measurements  along  a  transect  through  a  network  of  such  features.  A  probabilistic  determination  of  this 
error  is  described  here,  providing  a  starting  point  for  the  application  of  stochastic  geometry  theorems  in 
the  analysis  of  lead  geometries.  Errors  in  statistics  derived  for  other  lead  and  keel  features  are  discussed 
briefly.  While  the  application  is  to  sea  ice  leads  and  sonar  data,  the  methods  also  apply  to  the  general 
problem  of  sampling  linear  features  along  a  transect. 

1.4.1  DEFINITIONS,  NOTATION,  AND  AN  ILLUSTRATION 

In  the  following  discussions,  notation  follows  that  used  in  probability  theory,  where  Fz(:)  denotes 
the  distribution  function  (df)  for  the  population  random  variable  Z  with  specific  instance  z  (i.e.,  Fz(z)  s 
P[Z  <  z])  and  /z(z)  is  the  probability  density  function  (ndf).  Additionally,  E[Z]  and  Vc,.-[Z]  are  the 
expected  value  and  variance  of  Z. 

The  problem  is  to  relate  a  distribution  of  lead  widths  taken  along  a  line  perpendicular  to  the  local 
orientation  of  a  lead  (the  "actual"  width)  to  the  lead  widths  measured  along  a  transect  (the  "apparent"  lead 
width),  taking  into  account  lead  orientations,  and  lead  crossing  angles.  As  illustrated  in  Figure  1.27,  the 
following  continuous  random  variables  are  defined: 

X  -  actual  lead  width, 

X'  -  apparent  lead  width  measured  along  a  transect,  X'  S  X, 

0  -  lead  orientation  (0  2  0  <  7t), 

A  -  lead  intersection  angle  (0  £  A  <  n) 

with  specific  realizations  x,  xf,  0,  and  a.  Additionally,  let  <j>  be  the  transect  orientation  (0  £  <f>  <  rt).  The 
position  and  orientation  of  a  lead  within  the  plane  are  uniquely  specified  by  the  length  of  the  perpendicular 
that  connects  the  lead  to  the  origin,  and  the  angle  that  it  makes  with  a  fixed  reference  line.  The 
intersection  angle  A  is  measured  between  the  transect  and  the  lead,  anticlockwise,  and  is  the  difference 
between  their  orientations.  Finally,  define  A'  =  |n/2  -  A  |  to  be  the  crossing  angle  (ex'  in  Figure  1.27) 
measured  between  the  transect  and  a  perpendicular  to  the  lead  orientation  (0  £  A'  2  n/2).  The  relationship 
between  apparent  and  actual  lead  widths  is 


X  = 


X 

cos  (A') 


(12) 


where  X  2  X'.  Rearranging  terms,  a  lead  crossing  angle  can  be  determined  from  the  lead  widths  by 


A' 


=  COS'1 


The  potential  inaccuracies  of  measuring  lead  widths  along  a  transect  can  be  illustrated  by  randomly 
choosing  a  transect  orientation  and  location  on  a  satellite  image.  Here  we  provide  an  example  with  a 
Landsat  Multispectral  Scanner  (MSS)  band  4  (0.5-0.6  pm)  scene  of  the  Beaufort  Sea,  March  1988  (Figure 
1.28).  The  pixel  size  is  80  m;  image  size  is  80  x  80  km,  a  subset  of  a  Landsat  scene.  To  increase  the 
sample  size  of  lead  widths  measured  along  the  transect,  multiple  transects  of  the  same  orientation  are 
placed  randomly  on  the  image.  It  is  assumed  that  the  pattern  of  leads  is  similar  beyond  the  image 
boundaries.  Processing  of  the  Landsat  data  for  the  retrieval  of  lead  statistics  is  as  follows.  A  dynamic 
threshold  procedure  is  applied  that  estimates  the  probability  density  function  of  a  mixture  population 
(lead/ice)  for  small  regions  within  the  image,  and  a  binary  image  results.  Valid  lead  fragments  are 
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identified,  where  "valid"  refers  to  a  linear  feature  for  which  a  meaningful  width  and  orientation  can  be 
determined.  Linearity  is  determined  through  correlation/regression  analysis.  Lead  widths  are  measured 
perpendicular  to  the  regression  line,  every  kilometer  along  the  lead  length,  and  the  slope  of  the  regression 
line  is  a  measure  of  the  lead  orientation.  Further  details  of  this  procedure  are  given  in  Key  et  al.  [1990). 


Fig.  1.27.  The  geometrical  relationships  between  a  lead  and  a  transect.  See  text  for  definition  of  angles  and  length 
variables. 
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Fig.  1.28.  Landsat  MSS  band  4  scene  of  the  ice  pack  north  of  Alaska  in  March  1988.  Area  co/ered  is  approximately 
(80  km)2.  Field-of-view  is  80m. 


The  distribution  of  these  actual  lead  widths  x  is  shown  in  Figure  1.29  and  orientations  0  in 
Figure  1.30.  The  mean  lead  orientation  is  0.67  radians  (38°,  approximately  southeast  to  northwest  where 
the  top  of  the  image  is  north).  For  a  transect  orientation  «J>=3.0  radians  (172°,  south-southwest  to  north- 
northeast),  the  distribution  of  apparent  lead  widths  x!  is  illustrated  in  Figure  1.31,  with  crossing  angles  a' 
shown  in  Figure  1.32.  The  mean  actual  width  is  348  m  with  a  standard  deviation  of  201  m,  while  the 
mean  apparent  width  is  368  m  with  a  standard  deviation  of  474  m.  Additionally,  the  maximum  actual  lead 
width  in  the  image  is  1,376  m,  while  the  maximum  width  measured  along  the  transect  is  2,818  m.  With 
a  transect  orientation  of  0. 13  radians  (7.4°)  the  difference  between  the  actual  and  apparent  mean  widths 
is  139  m  and  the  maximum  width  is  2,670  m.  From  this  example  it  is  clear  that  significant  errors  can 
result  from  sampling  along  a  transect.  The  following  section  presents  a  method  to  assess  this  error. 
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1.4.2  PROBABILITY  MODELS 

Theorems  of  stochastic  geometry  that  are  applicable  uere  have  been  developed  through  the  study 
of  fibers  as  a  stationary  random  process  in  the  plane.  If  we  use  this  analogy  with  lead  networks,  then  after 
Stoyan  et  al.  (1987,  p.  240)  the  df  of  intersection  angles  is 

[*■“  sin  (6'  -  «  dFe( 6') 

W  -  4 -  (13) 

j,  [sin  (0'  -  «|  dFJQ') 

where  FA(a)  is  the  probability  of  intersection  angles  between  0  and  a,  Fe(0)  is  distribution  function  for 
lead  orientations,  dF0(Q)  =  /e(0)  aQ,  a  increases  in  an  anticlockwise  sense,  and  in  the  integral 
Fe(rc  +  a)  =  1  +  F0(a).  The  pdf /e(0)  may  be  an  assumed  mathematical  distribution  or  may  be  based  on 
an  observed  rose  of  directioa 

If  the  leads  are  isotropic  then  the  corresponding  orientations  have  a  simple  uniform  probability 
distribution  in  the  interval  0  £  0  <  n;  i.e.,  j0(Q)=l/n  for  all  0.  In  this  case  the  distribution  of  intersection 
angles  is  independent  of  the  transect  orientation.  The  probability  of  crossing  a  lead  that  is  oriented  across 
the  transect  (a  —>  n/2)  is  greater  than  for  one  running  more  parallel  (a  ->  0  or  re).  The  associated 
intersection  angles  have  density 

fA(a)  =  Vfc  sin  a  ,  O^a^n 

which  is  not  uniform  but  is  symmetrical  about  a  =  n/2.  The  corresponding  distribution  function  is 
Fa(«)  *  W)  -  (1  -  cos  a)  ,  O^oc^n 

which  is  a  special  case  of  (13)  for  /H(0)  =  l/n.  In  the  anisotropic  (preferred  orientation)  case  we  use  (13) 
for  the  distribution  of  intersection  angles,  and  the  corresponding  densities  are  determined  numerically. 

Two  different  intersection  angles  correspond  to  each  crossing  angle  so  that  the  distribution  and 
density  functions  for  the  crossing  angle  are 

FM)  s  W  * 

=  P[\l-A\  <.  a'] 

=  F[-a'  Z  (A-l)  ±  a'] 

=  F[(l-a')  <,  A  <,  (J4c0] 

=  FA(l+a')  -  FA(l-a') 

and 

-  ~  -  cO  * 

which  in  the  isotropic  case  yields  FA.(a')  =  sin  a'  and  /A.(a')  =  cos  a'. 
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An  expression  for  the  joint  density  function  of  the  apparent  and  actual  lead  widths  can  be  derived. 
Suppose  that  the  joint  pdf  of  X  and  A'  are  known,  which  will  b efXA(x,a')  =  fx(x)fA.(a')  if  the  two  variables 
are  independent.  Then  if  Y  and  Z  are  two  new  random  variables  that  are  functions  of  X  and  A'  such  that 
Y  =  X  and  Z  =  X'  =  X/cos  A',  then  the  joint  pdf  of  Y  and  Z  can  be  computed  using  a  standard  theorem 
(Ross,  1984,  p.  217): 


ty/y.4  -  >,Ax-X)  -  U*- cos  '(p)]  y 

=  ycos  '(y)]  y  [(*r  -  (14) 

The  first  of  these  expressions  is  valid  whether  or  not  X  and  A'  are  independent.  If,  however,  future 
research  indicates  that  large  leads  are  oriented  differendy  than  small  leads,  for  example,  then  the  joint 
density  function  must  be  determined  in  another  manner  (possibly  from  observations).  Using  (14)  the  pdf 
of  apparent  lead  widths  can  be  obtained: 


tJX)  *  i-Axy)  dx 

-  I'W  t'JOOS  '(y)]  y  -  X^)  dx  (,5) 

The  df  of  apparent  lead  widths  can  be  obtained  by  integrating  (15)  or  by  conditioning  on  the  value  of  X, 
again  assuming  that  X  and  A'  are  independent.  The  latter  method  yields 

FA*)  -  l'  FJcos-’f  *)]  fjx)  dx  (16) 

which  is  based  on  the  df  rather  than  the  pdf  of  A'.  The  functions  can  be  discretized  as  arrays  and  the 
integral  in  (16)  approximated  as  a  sum: 

Fr(y)  =  A  Y,  FM  W  •  '■  x^iA,  y-=jA  (17) 


where  FA.(j,i )  =  F A.[co$\ilj)\,  N  is  the  number  of  discrete  observations  and  A  is  the  increment  between 
observations.  If  these  functions  are  expressed  as  matrices,  (17)  becomes 


whose  derivation  is  given  in  Key  and  Peckham  (1991). 

In  this  study  the  error  in  measured  lead  width  is  defined  as  X'-X  (which  is  always  positive), 
although  other  definitions  such  as  X/X'  would  also  be  useful.  Equation  (14)  allows  us  to  compute  the 
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distribution  of  the  error  as  follows: 

FrJa)  =  P[X-X  £  a]  -  jfJT"  f„(x,X)  dx  dx ,  a  a  0 

For  the  isotropic  case 


M  =  I”  4(X)  sfa 


\j2x  +  a 


x  +  a 


dx 


All  moments  of  X'-X  can  be  computed  from/x..; 


Lead  Widths  in  Landsat  Image 


Width  (m) 


Fig.  1.29.  Lead  width  distribution  for  the  scene  in  Figure  1.28.  Widths  are  measured  along  a  perpendicular  to  the 
local  orientation  of  the  lead,  and  are  grouped  in  100  m  bins.  The  mean  width  is 
348  m,  the  standard  deviation  201  m,  and  maximum  width  1376  m. 
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Lead  Orientations  in  Landsat  Image 


Fig.  1.30.  Lead  orientations  for  the  scence  in  Figure  1.28.  The  mean  orientation  is  0.67  (38.4°)  with  standard 
deviation  0.87  (49.8°). 


Fig.  1.31.  Lead  widths  from  a  randomly-chosen  transect  across  Figure  1.28.  Transect  orientation  is  3.0'  (172°)  or 
approximately  south-southwest  to  north-northeast  where  the  top  of  tie  image  is  noth.  The  mean  width  is  368  m. 
the  standard  deviation  474  m,  and  maximum  width  2818  m. 
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Angle  (radians) 


Fig.  1.32.  Lead  crossing  angles  for  a  transect  across  Figure  1.28.  Transect  orientation  is  the  same  as  in  Figure  1.31. 
The  mean  crossing  angle  is  1.091  (62°)  with  standard  deviation  0.23'  (13°). 


1.4.3  APPLICATION 

These  models  are  now  applied.  First,  A  lead  width  distribution  measured  from  sonar  data  is  used 
to  estimate  the  actual  lead  width  distribution,  for  both  isotropic  and  anisotropic  orientations.  Lead 
orientation  and  actual  width  distributions  are  then  assumed  known,  and  the  expected  error  in  lead  width 
is  determined  for  a  variety  of  situations. 

Lead  width  distributions  have  been  described  by  power  laws  ( Wadhams ,  1981;  Steffen,  1987)  as 
have  floe  sizes  ( Rothrock  and  Thorndike,  1984).  The  negative  exponential  distribution  has  also  been  used 
( Dickins  et  al„  1986)  with  mean  lead  width  X  and  variance  X1.  The  exponential  model  implies  that  there 
are  a  finite  number  of  small  leads,  and  that  the  field  is  characterized  by  a  length  scale  X.  In  fact,  the  lead 
width  distribution  may  be  scale-free,  in  which  case  a  power  law  would  be  appropriate.  There  is,  of  course, 
a  Iowa  limit  imposed  by  the  resolution  of  the  measuring  instrument,  and  for  this  reason  as  well  as  for 
clarity  of  illustrating  expected  values,  we  use  the  negative  exponential  model 

Lead  orientations  may  be  random  or  may  have  a  preferred  orientation.  A  Gaussian  model  is  used 
here  for  preferred  orientations.  It  is  recognized,  howeva,  that  the  actual  shape  of  the  distribution  may 
be  bimodal,  whae  large  leads  with  one  orientation  are  intasected  by  smaller  leads  at  another.  Intosection 
angles  of  approximately  28°  have  been  observed  elsewhae  ( Marko  and  Thompson,  1977).  This  situation 
is  not  obvious  in  Figure  1.30,  although  the  distribution  is  not  purely  Gaussian  eitha. 

Table  1.6  lists  the  expected  error  for  a  variety  of  conditions,  whae  error  is  defined  by  the 
difference  between  the  actual  and  measured  lead  widths.  Case  1  considers  the  situation  where  the  apparent 
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lead  width  distribution  is  known.  The  apparent  lead  widths  are  based  on  submarine  sonar  data  recorded 
by  the  USS  QUEENFISH  in  August  of  1970  in  the  central  Canada  Basin  ( McLaren ,  1989).  Ice  draft  data 
were  measured  by  an  upward-beamed  fathometer  with  a  footprint  diameter  of  approximately  2.7  m  and 
a  vertical  accuracy  of  ±10  cm.  Sequences  of  continuous  points  with  drafts  230  cm  constitute  leads,  an 
example  of  which  is  given  in  Figure  1.33.  Given  a  mean  apparent  lead  width  of  60.6  m,  the  expected 
value  of  the  error  is  65.7  m  with  a  variance  of  1466.1  m1 2.  For  cases  2-7  in  Table  1.6  the  actual  width 
density  function  is  assumed  known,  and  the  apparent  lead  width  distribution  is  estimated.  For  cases  2  and 
3,  the  crossing  angle  distributions  are  shown  in  Figure  1.34,  and  the  error  distributions  in  Figure  1.35.  In 
the  preferred  orientation  cases  (2-5)  the  error  means  and  variances  are  clearly  dependent  upon  transect 
orientation. 

Other  applications  of  this  procedure  are  possible.  For  example,  laser  profilometer  transects  are 
analogous  to  sonar  transects,  and  the  methods  outlined  above  could  be  used  for  lead  and  ridge  spacing 
distributions  and  their  associated  errors.  As  in  the  illustration  with  Landsat  data,  transect  sampling  of 
satellite  imagery  is  a  natural  application.  Similarly,  heat  flux  through  leads  is  in  part  a  function  of  fetch, 
and  fetch  is  a  function  of  the  actual  lead  width  and  the  crossing  angle  of  the  wind.  If  the  wind  direction 
is  constant  as  it  travels  across  the  network  of  leads,  then  the  distribution  of  fetches  can  be  determined  from 
the  distribution  of  actual  lead  widths.  Finally,  it  may  be  possible  to  estimate  open  water  fraction  over  a 
large  area  from  the  apparent  lead  width  and  spacing  distributions  measured  along  a  transect.  This  research 
is  currently  in  progress,  with  results  to  be  presented  subsequently. 


Table  1.6.  Expected  error  in  lead  widths  (m)  under  a  variety  of  assumed  distributions  and  mean  values. 


Case 

fx 

fx- 

fs  ¥ 

E(X 

-X') 

Var(X  -  X') 

1 

?2 

Sonar 

Uniform 

•  a. 

65.7 

1466.1 

2 

A=20m 

? 

Gaussian3 

2.36  (135°) 

3.7 

3.1 

3 

A=20m 

? 

" 

0.52  (30°) 

32.5 

165.9 

4 

fc=40m 

? 

tt 

2.36  (135°) 

4.8 

6.4 

5 

A^40m 

? 

" 

0.52  (30°) 

36.0 

345.3 

6 

A^20m 

? 

Uniform 

— 

43.2 

653.1 

7 

A=40m 

? 

it 

” 

64.2 

1391.2 

1  Width  distribution  model  is  negative  exponential. 

2  "?"  refers  to  the  unknown  distribution. 

3  Parameters  of  the  Gaussian  model  are  p=n/4r  (45°)  and  a=0.3r  (17°). 
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Relative  Distance  (km) 


Fig.  1.33.  Submarine  sonar  ice  draft  data  for  a  2.5  km  section  within  the  Canada  Basin  north  of  Alaska.  Leads  are 
defined  as  continuous  sequences  of  points  with  drafts  no  greater  than  0.3  m  (dashed  line);  six  leads  occur  in  this 
section. 


Crossing  Angle  Distributions 

Transect  Orientations:  2.36  (solid),  0.52  (dashed) 


Crossing  Angle  (radians) 


Fig.  1.34.  The  distribution  functions  of  crossing  angles  FA.  for  cases  2  (solid)  and  3  (dashed)  in  Table  1.6. 


Probability 
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Error  Distributions 

Transect  Orientations:  2.36  (solid).  0.52  (dashed) 


Error  (m) 


Fig.  1.35.  The  distribution  functions  of  the  lead  width  errors.  Fx..x,  for  cases  2  (solid)  and  3  (dashed)  in  Table  1.6. 
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In  the  previous  section  sea  ice  leads  were  modeled  as  a  Poisson  line  process  and  a  method  was 
was  presented  for  estimating  the  "actual"  lead  width  distribution  from  a  width  distribution  measured  along 
a  transect  (the  "apparent"  lead  width  distribution).  The  purpose  of  this  section  is  to  extend  that  work  to 
the  estimation  of  the  fractional  area  coverage  of  leads  from  measurements  along  a  line.  This  line  could 
be  a  submarine  transect  under  the  ice  or  on  a  satellite  image.  Two  methods  from  stochastic  geometry  are 
presented  although  applications  have  not  yet  been  performed.  While  the  problem  is  illustrated  with  respect 
to  sea  ice  leads  it  should  be  noted  that  the  general  lineal  method  described  here  is  applicable  to  any 
geophysical  parameter  with  known  (or  determinable)  spatial  structure. 

1.5.1  GENERAL  LINEAL  METHOD 

The  general  expression  for  the  estimate  of  the  fractional  area  coverage,  /?',  of  a  geophysical 
parameter  whose  actual  fractional  coverage  is  p,  regardless  of  the  the  spatial  structure  of  that  parameter 
is 


pf  =  1 ±u  jj(x)dx. 


where  /(x)  is  the  indicator  function  for  the  underlying  function  q(\)  at  location  x  and  the  factor  is 
needed  for  normalization  and  depends  on  the  shape  of  the  structuring  element  U  (e.g.,  it  may  be  the  length 
of  a  line  or  the  area  of  a  square).  The  indicator  function  takes  on  a  value  of  1  if  q(x)  satisfies  some 
condition  and  0  otherwise.  For  example,  if  a  thresholding  procedure  is  used  to  determine  whether  or  not 
each  pixel  in  an  image  or  each  data  point  in  submarine  sonar  data  represents  some  phenomenon,  then 
/(x)=l  if  the  data  value  passes  the  threshold  test  and  /(x)=0  otherwise.  Again,  this  applies  to  any 
structuring  element  U. 

Following  Stoyan  et  al.  (1989)  the  expected  value  of  p\  E (p')>  is  p  and  its  variance  is 

Var(pO  =  EO tf-p)2  =  (Xy  jujuW  fix  _y  II)  dx  dy 


where  k,  is  the  autocovariance  function  oi  the  indicator  function  l: 

m  -  E[/(x)  /(x+/)J  -E[/(x)]2 


The  displacement  or  lag  r  =  |x-y|  is  such  that  the  covariance  depends  only  on  the  distance  between  the 
two  points  and  not  on  direction.  The  assumptions  are  that  the  geophysical  field  q(x)  is  stationary  and 
isotropic. 

Now  we  consider  the  case  where  the  structuring  element  is  a  line.  For  measurements  along  an 
array  of  N  parallel  lines,  each  of  length  L,  the  unbiased  estimate  of  the  fractional  area  coverage  is 


P 


I 

NL 


(18) 


where  /  is  the  total  length  of  NL  where  /(x)  =  1.  Extending  the  work  of  Rothrock  and  Thorndike  (1984), 
the  estimation  variance  is 
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Var(pO  =  N-'L-2  ^kl\x-y\)dxdy 

(19a) 

=  2N-'L2 ^L-t)k^dr 

(19b) 

where  Sf  is  the  set  containing  the  N  lines  and  r  =  |x-y|,  the  distance  between  two  locations  x  and  y  on 
the  line. 

Exponential  covariance  is  a  reasonable  model  for  many  geophysical  parameters  and  is  used  here: 

m  =  P(1  (r,a>  0) 

(20) 

where  a  describes  the  dependence  of  the  covariance  on  the  distance  r.  The  parameter 
determined  from  observed  autocovariances  by  rewriting  (20)  in  linear  form  as 

a  can  be 

In  [k{t)]  =  ln[p(1  -p)]  -ar 

(21) 

and  estimating  a  from  the  data  through  a  least  squares  regression. 

In  the  case  of  exponential  covariance  the  estimation  variance  in  (19)  is 

<a  _  1 

2p(1  -P)(  1  -  ,  ) 

Var  (pi)  =  aL 

aNL 

(22a) 

2p(1  -p)(  1  -_L) 

otL 

(22b) 

aNL 


1.5.2  APPLICATION 

One  real  and  two  simulated  images  are  used  in  the  application  of  the  above  methods.  The  real 
image  is  a  Landsat  MSS  band  4  (0.5-0.6pm)  scene  of  the  Beaufort  Sea,  March  1988.  The  binary  image 
produced  by  applying  a  threshold  to  the  original  grey-scale  image  is  shown  in  Figure  I.36a.  The  pixel 
size  is  80  m;  image  size  is  24  x  24  km,  a  subset  of  a  full  Landsat  scene.  Next,  a  network  of  leads  is 
simulated  as  a  Poisson  line  process.  The  mean  spacing  between  lines  (leads)  is  3000  m  and  their 
orientations  are  random.  The  lines  are  assigned  thicknesses  (widths)  following  the  negative  exponential 
density  function: 


W  -  -ie-*1 

where  w  is  lead  width  and  A.  is  the  mean  width.  For  the  simulation  A.  =  200  m.  One  realization  of  the 
Poisson  line  process  is  shown  in  Figure  I.36b,  again  as  a  binary  image,  where  the  pixel  size  is  137.5  m. 
Lastly,  a  cloud  field  is  simulated  as  an  ensemble  of  disks  whose  diameters  are  appioximately  normally- 
distributed  (in  a  true  Gaussian  distribution  negative  diameters  would  be  possible)  and  whose  center 
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locations  follow  a  binomial  point  process.  This  model  is  appropriate  for  cumuliform  clouds  but  is 
obviously  not  applicable  to  stratiform  cloud  decks.  One  realization  is  shown  in  Figure  I.36c. 

As  stated  earlier,  the  expected  value  of  p'  is  p;  i.e.,  the  mean  of  the  sampling  distribution  of 
sample  proportions  is  the  same  as  the  population  or  true  mean.  The  variance  of  the  sample  proportions 
is  given  by  (19)  in  the  general  case  and  (22)  for  exponential  covariance.  How  well  does  this  theory 
compare  to  observations?  As  the  first  step,  the  "true"  autocovariance  functions  were  estimated  from  six 
random,  horizontal  transects  through  each  image.  The  a  coefficient  in  (21)  was  computed  for  each  of  the 
six  transects  and  then  averaged.  With  the  simulated  clouds  and  leads  the  horizontal  transects  should 
adequately  represent  the  two-dimensional  structure  since  the  patterns  are  isotropic.  This  is  less  true, 
however,  with  the  Landsat  image  where  a  preferred  orientation  is  apparent.  Table  1.7  gives  the  results  of 
the  least  squares  fit  of  the  exponential  autocovariance  function  to  the  observed  autocovariances  in  the 
images.  Listed  are  the  regression-estimated  variance  (the  antilog  of  the  y-intercept  in  (21)),  a,  and  the 
correlation  coefficient  Next  the  distributions  of  sample  fractional  area  estimates  and  their  first  two 
moments  (mean  and  variance)  were  determined  by  computing  p'  with  (18)  for  each  of  500  single,  random, 
horizontal  transects  through  each  image  (i.e.,  the  number  of  transects  used  to  calculate  p'  in  (18),  N,  is  1 
for  each  of  500  samples).  Distributions  of  p'  were  also  computed  for  sets  of  ten  such  parallel  transects 
giving  5000  transects  (N=  10  for  each  of  500  samples).  Figure  1.37  shows  how  the  distribution  of  the 
estimates  changes  as  a  function  of  the  number  of  transects  or,  effectively,  total  transect  length.  For  single 
transects  (solid  lines)  a  broad  range  of  p'  are  possible. 
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rf~'.  Fig.  1.36.  (a)  Binary  image  based  on  a  Landsat  MSS 
band  4  scene  of  the  ice  pack  north  of  Alaska  in 
March  1988.  Field-of-view  is  80  m;  area  covered  is 
approximately  (24  km)2,  (b)  A  simulated  lead 

s  network  modeled  as  a  Poisson  line  process  with  thick 

- lines.  Field-of-view  is  137.5  m;  image  size  is 

approximately  (42  km)2,  (c)  A  simulated  cloud  field 

-1 -  based  on  a  random  disk  model.  Field-of-view  is 

137.5  m;  image  size  is  approximately  (42  km)2. 
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Table  1.8  gives  the  true  area  fracaon  p,  the  means  and  variances  of  the  estimated  area  fraction  p' 
from  the  observed  distributions  in  Figure  1.37,  and  the  variance  computed  from  (22b),  where  L= 304  pixels 
(same  units  as  r  in  (21)).  The  true  area  fraction  is  the  proportion  of  pixels  in  the  binary  image  where 
I(x)=  1.  It  is  apparent  that  the  theoretical  variance  of  the  estimate  given  in  (22)  is  generally  applicable  for 
the  simulated  leads  and  clouds  in  Figures  I.36b  and  I.36c.  For  the  Landsat  data,  however,  the  variance 
of  the  estimate  for  the  500-set  simulation  is  three  times  as  large  as  that  computed  with  (22)  or  a  factor 
of  1.7  for  the  standard  deviation.  This  is  due  to  the  anisotropic  nature  of  the  lead  network  in  the  image 
and  the  large  variability  in  the  autocovariance  function  computed  for  individual  transects.  A  two- 
dimensional  autocovariance  function  and  a  modification  of  (22)  may  be  needed  in  such  cases. 


Table  1.7.  Regression-estimated  parameters  of  the  autocovariance 
function  for  the  images  in  Figure  1.36. 


Figure 

Pd-P)1 
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R 
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-0.99 

lb 

lc 

-0.92 

1  Estimated  with  (21). 


Table  1.8.  Actual  and  estimated  fraction  area  coverages  for  Figures  I.36a-c 
using  one  and  ten  transects. 


Figure 

P 

N  =  1x500 

E  (p')  Var(p') 

N  =  10x500 

E  (//)  Vai(p') 

la 

G.035 

0.035  2.07e-3 

0.036  2.35e-4 

(6.84e-4)' 

(6.84e-5) 

lb 

0.067 

0.067  7.56e-4 

0.067  7.62e-5 

(7.38e-4) 

(7.38e-5) 

lc 

0.215 

0.218  7.28e-3 

0.214  7.78e-4 

(8.93e-3) 

(8.93e-4) 

1  Values  ir  parentheses  are  computed  with  (22b). 


The  fact  that  the  distributions  of  area  fraction  estimates  tends  toward  Gaussian  as  N  increases 
suggests  a  method  for  hypothesis  testing  and  confidence  interval  estimates.  If  a  normal  distribution  is 
assumed  to  apply,  then  the  probability  that  a  particular  area  fraction  estimate  comes  from  a  population 
with  area  fraction  p  can  be  determined.  This  is  perhaps  most  useful  for  confidence  interval  estimates  of 
the  true  area  fraction,  defined  as 
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pi  -  4P/2)  sd  {pi)  to  pf  +  z(P/2)  sd(p') 

where  sd(//)  =  [var(p')]w  is  the  standard  deviation  of  p'  and  1-p  is  the  level  of  confidence.  Since  the 
variance  depends  on  an  unknown  p,  p'  is  used  as  an  estimate  in  (22).  As  an  example,  suppose  that  for 
Figure  I.36b  samples  of  one  and  ten  transects  are  taken  and  for  both  p'= 0.05.  Assuming  the 
autocovariance  structure  given  in  Table  1.8  and  computing  the  variance  of  the  estimate  with  (22b)  the 
confidence  interval  estimate  of  the  true  area  fraction  at  the  90%  confidence  level  is  [0.011,  0.089]  for  a 
single  transect  and  [0.038, 0.062]  for  a  sample  of  ten  transects.  Neither  of  these  intervals  contains  the  true 
area  fraction  p=0.061.  Of  course,  the  probability  of  obtaining  such  a  p'  from  a  population  with  a  true 
fraction  of  0.067  is  very  small,  particularly  for  the  set  of  ten  transects  (0.01  as  opposed  to  0.24  for  a  single 
transect),  so  that  this  example  is  improbable  but  useful  for  illustration.  If,  on  the  other  hand,  we  obtain 
a  sample  p'  of  0.07  then  the  confidence  interval  estimate  of  p  is  [0.024, 0. 1 16]  for  a  single  sample  transect 
and  [0.056,  0.084]  for  a  set  of  ten.  Both  contain  the  true  area  fraction  but  the  larger  sample  size  gives 
a  much  smaller  range  for  the  estimate. 

The  shortcoming  of  this  approach  is  that  the  autocovariance  function  must  be  known.  It  may  be 
possible  to  estimate  it  from  the  data  itself  if  the  sample  size  is  large  enough,  although  this  is  somewhat 
circular.  In  some  cases  it  may  be  possible  to  infer  a  covariance  structure  by  assuming  a  simple  model  of 
the  geophysical  variable,  as  done  in  Rothrock  and  Thorndike  (1984)  for  sea  ice  floes.  Even  so,  some 
knowledge  of  the  field  is  needed;  in  their  case  the  diameter  of  the  floes.  If,  however,  some  basic 
autocovariance  structure  can  be  assumed  for  different  cloud  types,  sea  ice  leads,  etc.,  then  the  above 
procedure  is  certainly  useful  for  planning  sampling  studies,  id  probably  applicable  to  data  analysis  as 
well. 

1.5.3  SPECIAL  CASE;  POISSON  PROCESSES 

For  certain  stochastic  processes  it  is  possible  to  determine  the  fractional  area  coverage  from 
measurements  along  a  line  without  any  a  priori  knowledge  of  the  process.  Here  such  a  possibility  is  given 
for  a  Poisson  line  process  like  the  one  used  above  as  a  model  of  leads. 

For  a  Poisson  process  the  area  fraction  is  related  to  the  intensity1  t  of  the  process  and  the  mean 
"area"  of  the  objects 

pi  =  Prob[0  is  covered]  =  1  -e^  (23) 

where  O  is  an  arbitrary  origin.  The  area  measure  corresponds  in  units  to  the  intensity  me  e.g.,  for 
leads  the  intensity  is  the  number  of  points  per  unit  distance  and  £  is  mean  lead  width. 

The  area  fraction  can  be  now  estimated  from  lineal  measurements  through  the  use  of  the  line 
(lead)  thickness  (width)  distribution.  The  area  term  in  (23)  is  the  overall  mean  line  thickness,  VV,  defined 
as 


‘The  intensity  of  a  stochastic  process  is  commonly  called  the  density  of  the  process.  The  former  term 
is  used  here  in  order  to  avoid  confusion  with  the  concept  of  probability  density. 
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w  =  tt1  pi^de 

where  h<0)  is  the  mean  thickness  of  lines  with  orientation  9  (0  £  0  £  n)  ( Miles ,  1964).  This  applies  to 
lines  oriented  isotropically;  i.e.,  with  a  uniform  distribution  such  that /e(0)  =  n'1  where  f9  is  the  probability 
density  function  for  line  (lead)  orientations.  For  anisotropic  thick  lines  then 

w  -  £me)dF9(0) 


where  dF0(Q)  =/e(0)  dQ,  and  F0  is  the  cummulative  distribution  function  for  orientations.  A  method  for 
determining  the  actual  lead  width  distribution,  and  hence  W,  from  the  width  distribution  measured  along 
a  transect  has  been  presented  in  (Key  and  Peckham,  1991). 

As  an  example  of  the  use  of  (23),  the  lead  network  in  Figure  I.36b  was  generated  with  x  =  1/3 
(3  km  mean  spacing)  and  W  =  0.2  km.  This  gives  a  p'  estimate  using  (23)  of  0.064  compared  to  the  value 
of  0.067  reported  in  Table  1.8.  The  discrepency  is  a  function  of  the  image  creation  and  thresholding 
process,  where  all  leads  must  fill  an  entire  pixel. 

In  practice  the  intensity  of  the  process  is  not  known.  For  leads  modeled  as  a  Poisson  line  process 
an  estimate  of  x  can  be  obtained  from  the  transect  data,  where  the  points  of  intersection  of  the  transect 
with  leads  constitute  a  Poisson  process  of  intensity  2 z/n.  The  accuracy  of  this  estimate  depends  on  the 
size  of  the  region  over  which  the  measurements  are  made.  For  Figure  1. 36b  estimates  of  x  range  from 
0.19  to  upwards  of  0.45  which  results  in  an  estimate  of  p'  in  the  range  of  0.037  to  0.086.  There  is,  of 
course,  some  variability  in  the  estimate  of  W  as  well,  which  is  discussed  in  (Key  and  Peckham,  1991). 
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Part  ii:  Radiative  Transfer  Modeling  Studies 


II.l  RADIATIVE  TRANSFER  MODELING  AND  MODEL  VALIDATION 

In  order  to  utilize  radiative  transfer  models  most  effectively,  we  have  reviewed  how  they  treat 
factors  such  as  ice  cloud  morphology,  cloud  optical  thickness,  low  level  inversions,  boundary  layer  effects, 
and  other  aspects  relevant  to  the  polar  regions.  A  particular  concern  is  that  existing  models,  cloud 
properties,  standard  atmospheres,  etc.  have  typically  been  developed  for  low  and  mid-latitude  applications, 
and  may  thus  contain  biases  or  shortcomings  when  applied  to  polar  regions.  Part  of  this  effort  involves 
incorporating  polar  atmospheres  and  cloud  properties  into  the  models.  Arctic-specific  temperature  and 
humidity  profiles  have  been  obtained  and  are  discussed  in  the  next  section.  Unfortunately,  little 
information  on  the  microphysical  characteristics  of  arctic  clouds  is  available.  Data  that  are  available, 
however,  were  incorporated  into  the  models  where  appropriate.  This  is  limited  primarily  to  arctic  stratus 
experiments  during  the  early  1980’s,  and  some  measurements  of  aerosols.  Additional  data  is  expected  to 
be  obtained  during  LEADEX  by  instruments  on  the  NOAA  P-3. 

II.  1.1  RADIATIVE  TRANSFER  MODEL 

Work  described  in  this  part  of  the  report  relies  heavily  on  simulating  radiances  measured  by  the 
AVHRR  sensor.  To  simulate  radiances  in  the  AVHRR  thermal  channels,  daily  temperature  and  humidity 
profiles  in  each  season  are  used  with  the  LOWTRAN  7  radiative  transfer  model  {Kneizys  et  al.,  1988). 
Radiances  are  modeled  for  sensor  scan  angles  from  0°  to  60°  in  10°  increments.  The  appropriate  sensor 
response  function  is  applied  to  the  calculated  radiances,  and  radiances  are  then  converted  to  brightness 
temperatures.  Atmosperic  chemical  composition,  background  tropospheric  and  stratospheric  aerosols  for 
the  subarctic  winter  and  summer  models  are  used,  since  no  such  information  is  available  from  the  ice 
islands.  The  optical  properties  of  Arctic  haze  have  not  been  extensively  measured;  model  calculations 
(Blanche!  and  List,  1983)  show  that  the  volume  extinction  coefficient  of  Arctic  haze  is  generally  of  the 
same  order  of  magnitude  as  that  of  the  tropospheric  aerosols.  Therefore,  the  use  of  tropospheric 
background  aerosols  is  appropriate. 

In  order  to  test  some  sense  of  the  validity  of  the  radiative  transfer  model,  downwelling  longwave 
irradiances  (fluxes)  computed  with  LOWTRAN  were  compared  to  measurements  at  South  Pole,  Greenland 
and  Denver.  Radiosonde  data  from  the  three  locations  were  used  for  the  calculations.  To  obtain 
irradiance,  E,  from  LOWTRAN  (which  outputs  radiances,  L)  the  calculation  was  done  for  four  angles  (8 
streams)  and  employed  the  weighting  function; 

E  =  71(0.3626838^  +  0.3137066 L,  +  0.222381 0L,  +  0.1012285Z-4) 

where  the  subscripts  of  L  refer  to  the  sensor  zenith  angles  79.430°,  58.296°,  37.187°,  and  16.201”, 
respectively.  Additionally,  a  single-angle  method  (52.5°)  was  tested  and  found  to  give  fairly  accurate 
results.  The  bandwidth  of  3.5  -  50  pm  used  in  the  calculations  corresponds  with  the  bandwidth  in  which 
a  pyrgeometer  measures.  The  LOWTRAN  irradiances  differed  from  the  measurements  by  -4.7%  to  +5.6%. 
Assuming  5%  accuracy  for  the  pyrgeometer  data  these  results  are  acceptable.  Results  for  clear  sky  are 
listed  in  Table  II.l. 
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Table  n.  1.  Comparison  between  modeled  and  observed  longwave  fluxes  under  clear  sky  conditions. 


Station  and  Date 

Measured 
(W  nr2) 

4  Angles 
(W  m'2) 

Error  (%) 

52.5° 

(W  rn2) 

Error  (%) 

Denver,  9-30-89 

334 

336 

+0.60 

320 

-4.19 

South  Pole,  12-28-86 

128 

124 

-3.13 

117 

-8.60 

South  Pole,  7-31-86 

71 

75 

+5.63 

71 

0.0 

South  Pole,  11-6-86 

107 

102 

-4.67 

96 

-10.28 

Calculations  were  also  done  for  cloudy  sky  conditions.  In  comparisons  of  model  results  and 
observations,  the  observed  cloud  fraction,  Ac,  must  be  considered.  Here  the  model  is  used  to  estimate 
clear  and  overcast  fluxes,  and  these  are  weighted  by  the  observed  cloud  fraction. 

A  remaining  problem  concerns  the  differences  between  the  microphysical  properties  of 
LOWTRAN’s  cloud  models  and  those  of  the  observed  clouds,  which  are  unknown.  Table  n.2  gives  the 
results  of  the  cloudy  sky  comparisons. 

Table  II.2.  Comparison  between  modeled  and  observed  longwave  fluxes  under  cloudy  conditions. 

Station  and  Date 

Cloud  Fraction,  Type 

Measured 
(W  m'2) 

4  Angles 
(W  m'2) 

Error  (%) 

52.5° 

(W  nr2) 

Error  (%) 

Greenland,  7-23-90; 

8/8,  stratocumulus 

309 

301 

-2.59 

301 

-2.59 

Greenland,  7-26-90 

3/8,  stratocumulus 

283 

282 

-0.35 

275 

-2.48 

Greenland,  7-1-90 

1/8,  cirrus 

225 

227 

+0.89 

217 

-3.56 

South  Pole,  7-16-86 

8/8,  cirrus 

157 

156 

-0.64 

156 

-0.64 

South  Pole,  11-26-86 

8/8,  cirrus 

126 

129 

+2.38 

122 

-3.17 

To  examine  the  effect  of  vertical  temperature  structure  on  upwelling  longwave  radiation,  radiances 
in  the  three  channels  were  estimated  using  arctic  mean  and  subarctic  standard  winter  and  summer  profiles 
(described  below)  with  identical  seasonal  surface  temperatures.  The  maximum  difference  in  radiances  was 
0.05  W  m 2  sr'1  indicating  that  the  vertical  temperature  distribution  of  the  relatively  dry  arctic  atmosphere 
plays  a  very  small  role  in  the  attenuation  of  upwelling  longwave  radiation. 
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II.  1.2  ICE  ISLAND  DATA 

Analyses  are  based  on  atmospheric  temperature  and  humidity  profiles  collected  by  rawinsonde 
from  a  Soviet  ice  island  (NP-26),  located  at  approximately  85°  N  170°  W  during  1983-1987  (Figure  II.  1). 
Generally  two  observations  per  day  were  collected  covering  a  vertical  range  of  0-25  km.  Profiles  that 
have  at  least  10  levels  are  retained  in  the  analysis.  Observations  include  temperature,  dew  point 
depression,  wind  speed,  and  wind  direction.  For  the  years  1986-87  surface-based  cloud  observations  are 
also  available.  These  observations  include  low,  middle,  and  high  cloud  types,  height  of  the  cloud  base, 
and  cloud  fraction. 

Only  clear  sky  profiles  are  of  interest  in  one  of  the  studies  below,  and  since  the  satellite  thermal 
radiances  under  cloudy  conditions  will  reflect  cloud  top  temperature  and  a  significant  amount  of  cloud 
cover  will  affect  the  lower  tropospheric  temperature  structure,  clear  sky  "seasons"  that  differ  in  their 
vertical  temperature  and  humidity  structures  are  then  defined.  The  seasons  are  determined  objectively  with 
a  squared  Euclidean  distance  clustering  algorithm;  the  variables  are  temperature  and  humidity  at  each 
level.  To  reduce  the  degree  of  statistical  dependence  between  levels,  only  one  measurement  per  kilometer 
was  used.  The  resulting  seasons  are  winter:  October  through  March,  summer:  June  through  August,  and 
transition:  April,  May,  and  September.  The  resulting  mean  seasonal  temperature  profiles  for  clear,  cloudy 
(greater  than  75%  cloud  cover),  and  mixed  conditions  are  shown  in  Figure  II.2. 


Figure  II.  1.  Location  of  the  Soviet  ice  island  NP-26. 


Temperature  (K) 


Figure  n.2.  P  otic  winter,  transition,  and  summer  temperature  profiles  under  clear,  cloudy,  and  mixed  conditions 
from  a  Soviet  ice  island  located  in  the  Canada  Basin. 
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Here  we  address  the  need  to  understand  how  different  sensors  respond  radiometrically  to  varying 
surface  types  and  intervening  atmospheric  layers  by  evaluating  the  combined  use  of  multispectral  image 
and  contrast  analyses  to  determine  thresholds  of  lead  detection  under  varying  atmospheric  conditions.  Our 
objective  is  to  estimate  the  narrowest  resolvable  lead  under  a  given  set  of  surface  and  atmospheric 
conditions,  sensor  FOV  and  viewing  geometry.  Although  this  investigation  focuses  on  the  use  of  AVHRR 
data,  the  approach  should  have  general  applicability  to  other  thermal  sensors  of  differing  spectral  response 
and  view  characteristics.  The  modeling  results  presented  here  contrast  simulations  of  dear-sky  conditions 
with  those  including  the  effects  of  various  types  of  horizontally  homogeneous  cloud  or  aerosol  layers. 
Cloud  detection  using  AVHRR  data  is  not  addressed.  For  a  discussion  of  polar  cloud  detection  see  Key 
and  Barry  (1989),  Sakellariou  et  al.  (1992)  and  Yamanouchi  et  al.  (1987)  and  references  therein. 

Our  overall  approach  is  as  follows.  We  first  simulate  top-of-the-atmosphere  (TOA)  radiances  for 
three  thermal  channels  of  the  AVHRR  instrument.  Average  clear-sky  January  conditions  for  the  central 
Arctic  are  assumed.  Operating  in  the  NIR  and  IR  "atmospheric  window"  regions  of  the  spectrum,  these 
channels  are  especially  sensitive  to  surface  emissions,  but  are  also  affected  by  any  intervening  atmospheric 
layer  that  absorbs/emits  thermal  radiation.  Thus,  surface  (or  skin)  temperatures  and  emissivities  are  varied 
to  evaluate  the  effects  of  different  surface  types.  All  modeled  radiances  are  convened  to  equivalent 
blackbody  temperatures  (brightness  temperatures)  by  inverting  the  Planck  function  ( NOAA ,  1991)  to 
facilitate  comparisons  with  physical  temperatures  and  the  analyses  of  bispectral  results  (NIR-IR  brightness 
temperature  differences)  and  derived  thermal  contrasts.  Next,  model  clouds  or  haze  layers  are 
hypothetically  inserted  into  the  atmosphere  to  examine  the  behavior  of  simulated  brightness  temperatures 
and  brightness  temperature  differences  (BTDs)  as  a  function  of  layer  optical  depth.  These  "split  window" 
results  are  then  examined  for  signatures  that  characterize  haze,  stratiform  water  clouds,  clear-sky  ice 
crystal  precipitation  (ICP)  or  high  level  cirrus  clouds.  Finally,  by  normalizing  the  difference  between 
channel  brightness  temperatures  of  a  lead  pixel  and  its  background  (i.e.  the  multiyear  ice  pack)  by  the 
brightness  temperature  of  the  background  scene  normalized  contrast  values  are  derived  and  evaluated  as 
a  means  to  determine  the  limits  of  lead  detection  given  certain  sensor  characteristics  and 
atmospheric/surface  properties. 

II.2.1  RADIANCE  SIMULATIONS 

The  radiative  transfer  code  LOWTRAN  7  ( Kneizys  et  al.,  1988)  (hereafter  simply  LOWTRAN) 
is  used  to  compute  top-of-the-atmosphere  upwelling  radiances  from  which  satellite-derived  brightness 
temperatures  are  simulated.  In  cases  involving  haze,  changes  in  the  properties  of  aerosols  as  a  function 
of  relative  humidity  (RH)  are  accounted  for  (e.g.,  Shettle  and  Fenn,  1979;  Blanchet  and  List,  1987)  by 
first  modifying  the  effective  refractive  indices  of  the  bimodal  particle  size  distribution  and  then 
recomputing  extinction  and  absorption  coefficients  based  on  Mie  theory  ( Kneizys  et  al.,  1980). 

In  this  study  we  simulate  AVHRR  radiances  for  channels  3, 4  and  5.  Channel  3  measures  in  the 
NIR  window  region  of  the  spectrum  and  is  centered  at  3.7  pm  while  the  IR  channels  4  and  5  are  centered 
at  10.8  pm  and  12.0  pm,  respectively.  All  our  calculations  include  the  effects  of  multiple  scattered 
thermal  radiation  and  are  made  at  steps  of  5  cm'1  (equivalent  to  0.06  pm  at  11  pm).  The  angular  and 
spectral  dependencies  of  snow  and  water  emissivities  are  also  taken  into  account.  In  all  cases, 
LOWTRAN  was  initialized  for  average  clear-sky  January  temperature  and  humidity  profiles  based  on  an 
analysis  of  Soviet  ice  island  data  collected  in  the  central  Arctic  (e.g..  Key  and  Haefliger,  1992;  Serreze 
et  al.,  1992).  The  mean  clear-sky  January  temperature  and  dewpoint  temperature  profiles  are  shown  in 
Figure  II.3  with  the  boundaries  of  subsequently  prescribed  hypothetical  layers  of  haze,  stratus  cloud,  cirrus 
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cloud  and  ice  crystal  precipitation  indicated.  Because  there  is  essentially  no  information  available  on  the 
vertical  structure  of  atmospheric  gases  in  the  central  Arctic,  we  assume  that  average  "subarctic  winter" 
concentrations  of  03,  CH4,  CO  and  NaO  exist  in  the  atmosphere  when  running  LOWTRAN.  Model 
subarctic  winter  background  aerosol  concentrations  for  the  troposphere  (2  to  10  km)  and  the  stratosphere 
(10  to  30  km)  were  also  prescribed.  To  simulate  haze  effects,  boundary  layer  (0  to  2  km)  aerosol 
concentrations  were  varied  by  specifying  the  layer  visibility,  but  in  all  other  cases  the  default  "rural" 
aerosol  model  for  the  boundary  layer  was  used.  To  examine  the  full  range  of  AVHRR  scan  angles  (0° 
to  approximately  55°)  we  made  calculations  at  0°  (nadir),  20°  and  50°.  Only  results  for  0°  and/or  50°  are 
presented  here. 


Mean  January  Temperature 
Clear  Sky 


Fig.  n.3.  Mean  clear-sky  January  temperature  and  dewpoint  temperature  profiles  for  the  central  Arctic.  Also  shown 
are  the  vertical  positions  of  hypothetical  layers  of  cirrus  cloud,  boundary  layer  haze,  ice  crystal  precipitation  and 
stratus  cloud  that  are  prescribed  for  radiative  transfer  simulations  that  are  described  in  the  text. 
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II.2.2  INTERVENING  ATMOSPHERIC  EFFECTS 

The  effects  of  four  commonly  observed  atmospheric  phenomena  in  the  Arctic  are  considered:  1) 
Arctic  haze  (e.g.,  Rahn  and  McCaffrey,  1980;  Valero  et  al.,  1989;  Valero  et  al.,  1983)  which  persists  over 
large  regions  of  the  Arctic,  especially  during  late  winter  and  early  spring,  2)  Arctic  stratus  cloud  (e.g., 
Tsay  and  Jayaweera,  1984;  Tsay  et  al.,  1989;  Curry  et  al.,  1992)  which  commonly  obscures  satellite 
viewing  of  the  sea  ice  surface  in  summer,  but  may  be  thin  enough  during  the  winter  months  to  make  lead 
detection  feasible,  3)  dear-sky  ice  crystal  precipitation  (ICP)  which  has  a  significant  effect  on  the  radiation 
balance  of  the  surface/atmosphere  system  in  the  Arctic  ( Curry  et  al.,  1990;  Curry  et  al.,  1989a;  Curry  et 
al.,  1989b)  and  4)  high  level  cirrus  clouds  which  are  the  "most-frequently-occurring  doud  type"  observed 
in  the  central  Arctic  during  winter  and  spring  ( Warren  et  al.,  1988).  Each  of  these  has  a  distinct  effect 
on  the  upwelling  thermal  radiation  emitted  from  the  underlying  surface  and  atmosphere  depending  on  their 
microphysical  properties,  geometric  thicknesses  and  positions  within  the  atmosphere.  The  differences  in 
the  radiative  properties  of  atmospheric  aerosols  (haze),  water  droplets  or  ice  crystals  result  in  varying 
degrees  of  scattering  and  absorption  as  a  function  of  wavelength.  These  differences  can  be  exploited  using 
multispectral  techniques  to  distinguish  various  types  of  attenuating  layers  that  may  exist  in  the  Arctic 
atmosphere  assuming  that  the  underlying  surface  properties  can  be  determined  by  other  means. 

11.2.2.1  Arctic  Haze 

The  optical  properties  of  Arctic  haze  have  not  been  extensively  measured,  but  model  calculations 
indicate  that  the  volume  extinction  coefficients  of  Arctic  hazes  are  to  a  first  approximation  the  same  order 
of  magnitude  as  those  for  tropospheric  aerosols  (e.g.  Blanc het  and  List,  1983;  Tsay  et  al.,  1989).  Because 
Arctic  haze  generally  contains  an  anthropogenic  component  of  carbonaceous  material  transported  from  the 
lower  latitudes  ( Rosen  et  al.,  1981;  Kahl  and  Hansen,  1989),  the  "urban"  aerosol  model  of  LOWTRAN 
was  selected  to  simulate  low  level  haze  layers.  This  model  represents  a  mixture  of  20%  soot-like  aerosols 
and  80%  rural  type  aerosols  contained  in  the  0  to  2  km  boundary  layer  ( Kneizys  et  al.,  1980).  The 
extinction  coefficient  p  for  boundary  layer  haze  as  defined  in  LOWTRAN  is  determined  from  • .  prescribed 
atmospheric  visibility  V  using  Koschmieder’s  formula:  V  =  1/p  ln(l/6),  where  6  is  a  threshold  contrast 
taken  to  be  0.02. 

The  infrared  opacity  of  aerosol  layers  is  known  to  increase  quite  dramatically  with  increasing 
relative  humidity  ( Blanchet  and  List,  1987;  Shettle  and  Fenn,  1979),  thus  an  assessment  of  how  water 
uptake  by  hygroscopic  aerosols  affects  simulations  of  brightness  temperatures  and  BTDs  was  also  made. 
Results  for  a  saturated  haze  layer  (RH  =  99.9%)  composed  of  "wet"  aerosol  particles  are  contrasted  with 
those  for  moderately  dry  (RH  =  70%)  haze  layers  found  to  characterize  mean  January  conditions  in  the 
Arctic.  LOWTRAN  is  designed  to  modify  the  absorption  and  scattering  coefficients  of  aerosol 
distributions  by  1)  assuming  growth  of  particulates  as  a  function  of  RH  based  on  the  results  of  Hanel 
( 1976),  2)  adjusting  their  effective  refractive  indices  and  3)  recomputing  their  radiative  properties  based 
on  Mie  theory  ( Kneizys  et  al.,  1980). 

11.2.2.2  Arctic  Stratus 

Arctic  cloud  climatologies  show  marked  increases  in  average  low  cloud  amounts  during  spring 
attributed  to  the  presence  of  stratiform  water  clouds  within  the  boundary  layer  which  reach  a  maximum 
coverage  of  about  70%  during  the  summer  months  ( Huschke ,  1969;  Vowinckel  and  Orvig,  1970).  Stratus 
clouds  have  a  potentially  dramatic  impact  on  the  surface-atmosphere  heat  budget  depending  on  whether 
their  shortwave  albedo  effects  or  longwave  greenhouse  effects  dominate  the  lower  tropospheric  radiation 
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balance  ( Curry  et  al.,  1992).  Summertime  visible  stratus  optical  depths  tend  to  be  large,  in  cases 
exceeding  20  ( Herman  and  Curry ,  1984),  thus  the  use  of  visible  imagery  to  detect  leads  is  not  practical, 
and  because  thermal  contrasts  between  leads  and  multiyear  sea  ice  tend  toward  zero  as  cloud  opacity 
increases,  detection  using  thermal  imagery  is  also  impractical  in  most  cases.  During  winter,  however, 
stratus  cloud  cover  is  often  less  than  20%  and  the  cloud  layers  tend  to  be  thin  optically  with  visible  optical 
depths  on  the  order  of  2  ( Curry  et  al.,  1992).  Detecting  leads  using  thermal  imagery  in  winter  may, 
therefore,  be  possible.  Using  the  LOWTRAN  model  for  "stratus"  we  assess  the  feasibility  of  detecting 
leads  dining  winter  by  assuming  that  a  horizontally  homogeneous  cloud  layer  380  m  thick  exists  in  the 
lower  atmosphere  as  indicated  in  Figure  1.  For  stratus  simulations  the  desired  range  of  visible  (0.55pm) 
optical  depth  x  was  obtained  by  varying  the  conversion  factor  from  equivalent  LWC  (g  m 3)  to  extinction 
coefficient  (km'1)  in  the  LOWTRAN  code  ( Kneizys  et  al.,  1988).  The  stratus  droplet  size  distribution  is 
represented  by  a  modified  gamma  distribution: 

n(t)  =  27  r2e-°6r 

where  r  is  the  droplet  radius  and  n  is  number  density.  The  total  number  density  is  taken  to  be  250  cm'3 
and  the  mode  radius  is  6.67pm  (for  mass  distribution).  Details  of  the  LOWTRAN  cloud  models  are  given 
in  Shettle  et  al.  (1988). 

11.2.2.3  Layers  of  Ice  Crystals 

Both  the  high  level  cirrus  cloud  and  low  level  ice  crystal  precipitation  simulations  were  made  by 
inserting  the  LOWTRAN  "standard"  cirrus  model  (, Shettle  et  al.,  1988)  into  the  atmosphere  at  the  vertical 
positions  indicated  in  Figure  II.3.  For  both  conditions  the  desired  range  of  optical  depth  was  obtained  by 
assuming  the  appropriate  values  of  0.55pm  volume  extinction  coefficients  for  a  2  km  thick  cirrus  based 
at  8  km  and  a  1  km  thick  ICP  layer  based  at  the  surface.  Note  that  ice  crystal  precipitation  has  been 
observed  from  the  surface  to  heights  exceeding  3  km  in  the  Arctic,  but  it  is  most  frequently  observed 
below  about  l  km  (Curry  et  al.,  1990). 

For  all  of  the  cases  discussed  above  theoretical  calculations  were  made  for  a  visible  optical  depth 
range  of  0  (clear-sky)  to  100,  but  results  are  presented  only  for  values  between  0  and  10.  Realistic  layer 
thicknesses  and  extinction  values  of  aerosols  and  ice  crystals  are  such  that  visible  optical  depths  rarely 
exceed  10  and  although  stratus  optical  depths  may  exceed  this  value  during  the  summer  months,  low 
thermal  contrasts  and  high  visible  opacity  preclude  lead  detection  at  this  time  of  year.  In  the  winter,  even 
when  mixed-phase  layers  contain  small  amounts  of  liquid  water,  optical  depths  are  generally  within  the 
range  represented  here. 

11.2.3  SURFACE  CHARACTERISTICS 

Model  runs  were  initialized  for  three  different  surface  temperatures  to  characterize  open  or 
refrozen  leads  and  a  fourth  temperature  representing  the  surrounding  ice  pack  which  is  assumed  to  be  2 
m  thick  and  in  equilibrium  with  the  surface  air  temperature.  In  the  discussion  that  follows,  the  terms 
"skin"  and  "surface"  temperature  are  used  interchangeably  and  should  not  be  confused  with  shelter 
temperature  (measured  2  m  above  ground  level)  which  is  generally  higher  than  the  actual  skin  temperature 
when  a  surface-based  temperature  inversion  exists.  Shelter  temperature  may  differ  from  skin  temperature 
by  more  than  10°C  depending  on  the  region  and  time  of  year  (e.g.,  Stowe  et  al.,  1988,  Rossow  et  al., 
1989).  The  following  temperatures  were  computed  to  represent  a  range  of  lead  types:  for  open  leads,  27 1 
K;  for  leads  covered  by  5  cm  thick  ice,  256  K;  and  for  15  cm  thick  leads,  248  K.  The  ice  pack  was 
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assumed  to  be  in  thermal  equilibrium  with  the  surface  air  and  was  assigned  the  temperature  of  235  K 
corresponding  to  the  first  rawinsonde  level. 

All  surfaces  with  the  exception  of  open  leads  are  assumed  to  be  snow-covered  and  directional 
snow  emissivities  were  modeled  following  the  procedure  in  Dozier  and  Warren  (1982).  This  method 
involves  calculating  the  single  scattering  albedo,  asymmetry  factor  and  phase  function  for  snow  grains 
using  a  Mie  code  before  determining  the  directional,  wavelength-dependent  emissivities  using  the  delta- 
Eddington  approximation  of  the  radiative  transfer  equation.  These  emissivities  were  then  integrated  over 
the  response  function  for  each  channel  i: 

[Xl  e  (A.0) <|>,(*)  dk 

t  ((0)  =  A 


where  e(X,0)  is  the  emissivity  in  direction  0  at  wavelength  X  and  <j>j  is  the  ilh  sensor  response  functioa 
Key  and  Haefliger  (1992)  note  that  differences  between  the  integrated  emissivities  for  the  NOAA  series 
7,  9,  and  1 1  AVHRR  sensors  are  on  the  order  of  0.0001.  In  the  current  study  we  use  the  NOAA  7  values 
only.  For  brightness  temperature  simulations  over  open  leads,  the  angular  dependence  of  channel 
emissivities  of  water  are  determined  through  Fresnel  calculations.  Table  II.3  gives  the  angular  emissivities 
used  in  the  current  analysis  as  a  function  of  AVHRR  scan  angle  and  surface  type.  In  reality,  newly 
refrozen  'eads  are  clear  of  snow,  thus  pure  ice  emissivities  should  be  used  for  best  results.  However,  the 
authors  are  aware  of  no  comprehensive  set  of  measurements  nor  method  from  which  the  spectral, 
directional  emissivities  of  planar  ice  can  be  determined.  The  emissivity  of  ice  approaches  unity  and  is 
often  assumed  to  be  1.0  for  field  investigations  (e.g.,  Konig-Langlo  andZachek,  1991).  A  more  accurate 
value  may  be  0.97  ( Hobbs ,  1974)  which  is  within  2%  of  the  snow  values  listed  in  Table  II.3.  Thus,  the 
assumption  that  all  surfaces  are  snow-covered  will  not  result  in  serious  errors  in  simulated  brightness 
temperatures  nor  will  differences  in  channel  brightness  temperatures  or  contrast  ratios  be  affected 
significantly. 


TABLE  II.3.  Angular  emissivities  of  snow 
and  water  in  NOAA  7  AVHRR  channels  3, 
4  and  5  at  two  satellite  scan  angles. 


Snow 

Water 

Channel 

0° 

50° 

0° 

LA 

O 

o 

3 

.998 

.992 

.976 

.961 

4 

.999 

.996 

.992 

.984 

5 

.996 

.987 

.986 

.972 

For  even  more  accurate  simulations  of  brightness  temperatures,  channel  radiances  should  also  be 
integrated  over  the  appropriate  sensor  response  functions.  These  vary  from  one  satellite  to  another  as 
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discussed  in  Key  and  Haefliger  (1992).  In  that  study,  they  found  that  by  using  the  full  response  functions 
instead  of  assuming  rectangular  functions  (i.e„  100%  response  at  all  wavelengths  within  the  portion  of  the 
channel  where  the  actual  response  is  at  least  50%)  the  brightness  temperature  differences  were  only  on 
the  order  of  0.05  K  for  channel  4,  but  were  about  0.5  K  for  channel  5  assuming  typical  January  conditions 
in  the  central  Arctic.  Because  we  focus  here  on  channel  4  results,  we  use  the  rectangular  function  for  all 
simulations  rather  than  performing  this  time-consuming  integration.  This  is  justified  in  that  the 
computation  of  normalized  contrasts  are  based  on  differences  and  ratios  of  brightness  temperatures  at  one 
or  another  wavelength  so  that  small  absolute  biases  due  to  differences  in  response  functions  have  a 
negligible  effect  on  the  final  results. 


II.2.4  RESULTS 

II.2.4.1  Simulated  Brightness  Temperatures 

Examples  of  channel  brightness  temperatures  as  a  function  of  visible  (0.55  pm)  optical  depth  are 
shown  in  Figure  II.4  for  boundary  layer  haze,  ice  crystal  precipitation  and  high-level  cirrus  clouds  at  0° 
and  50°  satellite  scan  angles  for  channels  3,  4,  and  5  of  the  AVHRR  instrument.  These  layers  were 
positioned  as  shown  in  Figure  1  and  mean  January  temperature  and  humidity  profiles  were  assumed. 
Within  each  panel  (from  top  to  bottom)  are  plots  that  relate  to  the  different  prescribed  surface  temperatures 
representing  open  leads,  leads  of  5  cm  and  15  cm  thickness  (271  K,  256  K  and  248  K,  respectively),  and 
the  multiyear  pack  ice  (235  K).  Note  that  the  boundary  layer  is  moderately  dry  during  January  in  the 
central  Arctic  with  RH  averaging  about  70%  between  the  surface  and  2  km  suggesting  a  predominance 
of  clear  skies  containing  low  concentrations  of  dry  aerosols  during  this  time  of  the  year.  At  zero  optical 
depth  the  physical  surface  temperatures  are  reasonably  well  represented  by  simulated  channel  brightness 
temperatures  because  the  selected  channels  are  all  within  NIR  and  IR  window  regions  of  the  spectrum 
where  sensitivity  to  the  relatively  dry  Arctic  atmosphere  is  least  and  because  open  water,  newly  refrozen 
leads  and  snow-covered  surfaces  all  have  high  emissivities.  Regardless  of  the  underlying  surface  type, 
channel  3, 4  and  5  brightness  temperatures  tend  to  converge  to  the  blackbody  radiating  temperature  of  the 
top  of  the  intervening  cloud  or  haze  layers  as  optical  depths  increase,  though  the  rate  of  convergence 
varies  depending  on  the  microphysical  properties  of  the  intervening  layer  and  the  layer’s  position  and 
mean  temperature  relative  to  the  surface.  Scan  angle  effects  are  also  apparent  by  comparing  corresponding 
panels  for  0°  and  50°  viewing  angles.  No  matter  what  type  of  surface  or  intervening  medium  exists,  the 
convergence  of  simulated  brightness  temperatures  to  layer  top  temperatures  occurs  faster  when  viewing 
off-nadir.  This  is  due  to  the  increase  in  optical  path  length  by  a  factor  l/cos(0),  where  0  is  the  satellite 
scan  angle  as  viewing  angles  increase  from  nadir. 
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Fig.  n.4.  Simulated  brightness  temperatures  for  layers  of  boundary  layer  haze,  ice  crystal  precipitation  and  high-level 
cirrus  cloud  for  three  AVHRR  thermal  channels  assuming  four  surface  types  (from  top  to  bottom  at  left  of  each 
panel):  open  water,  5,  15,  and  200  cm  thick  ice.  Results  are  shown  for  satellite  view  angles  of  0“  and  5 O'  over  a 
range  of  0.55pm  optical  depths  between  0  and  10.  Mean  clear-sky  January  temperature  and  humidity  profiles  for 
the  central  Arctic  are  assumed. 
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Because  haze  particles  absorb  and  emit  much  less  NIR  and  IR  radiation  than  do  ice  panicles  or 
water  droplets  at  a  given  visible  optical  depth,  haze  layers  neither  attenuate  nor  enhance  significantly  the 
thermal  emissions  from  the  surface.  Thus,  simulated  NIR  and  IR  brightness  temperatures  during  typical 
winter  conditions  are  not  very  sensitive  to  changes  in  aerosol  loading  in  the  boundary  layer.  In  these 
situations,  the  relatively  dry  haze  layers  become  essentially  black  at  visible  optical  depths  exceeding  95 
(equivalent  to  an  infrared  optical  depth  of  about  6).  This  extreme  value  corresponds  to  an  unrealistic 
visibility  of  less  than  85  m.  As  noted  by  Blanchct  and  List  (1987),  however,  the  infrared  opacity  of 
aerosols  increases  dramatically  as  relative  humidity  increases.  In  fact,  IR  optical  depths  may  exceed  those 
for  visible  wavelengths  in  saturated  haze  layers.  This  phenomenon  occurs  due  to  the  uptake  of  available 
water  by  hygroscopic  particles  such  as  sulfuric  acid  and  deliquescent  compounds  within  the  layer  causing 
a  shift  to  a  larger  particle  size  distribution  and  corresponding  enhancements  of  the  absorption  and 
scattering  cross  sections  of  the  particles. 

We  wvaluated  the  effects  of  increasing  relative  humidity  on  aerosol  layers  by  hypothetically 
saturating  the  0  to  2  km  layer  and  recalculating  the  brightness  temperatures  for  a  "wet"  aerosol  layer.  This 
is  easily  accomplished  using  LOWTRAN  because  the  absorption  and  scattering  coefficients  are  modified 
in  accordance  with  hygroscopic  growth  as  described  previously.  The  optical  depth  dependence  of  channel 
4  brightness  temperatures  for  a  haze  layer  with  mean  RH  =  99.9%  was  found  to  be  virtually  identical  to 
that  for  a  haze  layer  having  RH  =  70%  while  channel  5  brightness  temperatures  tended  to  converge 
slightly  faster  to  the  blackbody  temperature  of  the  layer.  The  NER  channel  3  simulations  revealed  the  most 
pronounced  change  due  to  saturation.  They  indicated  slightly  faster  convergence  with  increasing  optical 
depth  than  either  of  the  IR  channels  and  also,  interestingly,  NIR  values  converge  to  a  temperature  below 
the  physical  temperature  of  the  layer  top.  This  leads  to  negative  NIR-IR  spectral  differences. 
Qualitatively,  an  optically  thick,  saturated  haze  layer  has  similar  radiative  behavior  as  does  a  stratiform 
water  cloud.  As  these  layers  increase  in  opacity  they  become  black  to  infrared  radiation  resulting  in  a 
rapid  convergence  of  brightness  temperatures  (no  matter  what  the  underlying  surface  is)  to  the  physical 
temperature  of  the  layer  top;  but  in  the  NIR,  thick  layers  totally  attenuate  the  upwelling  radiation  from 
the  surface  while  contributing  little  to  TOA  radiances  because  these  layers  have  low  NIR  emissivities. 
Thus,  NIR  brightness  temperatures  are  actually  colder  than  corresponding  IR  values  (e.g.,  Yamanouchi  et 
al.,  1987).  These  are  important  considerations  because  the  relative  magnitudes  of  NIR-IR  bispectral 
differences  as  a  function  of  relative  humidity,  phase  and  optical  depth  may  be  exploited  to  distinguish 
different  types  of  intervening  atmospheric  layers. 

The  importance  of  knowing  the  vertical  position  of  an  intervening  layer  is  apparent  by  comparing 
the  results  for  the  low  level  ICP  with  those  for  high  cirrus  clouds.  Both  layers  have  identical 
microphysical  and  radiative  properties  based  on  the  standard  LOWTRAN  cirrus  model.  However,  their 
net  radiative  effect  on  the  upwelling  radiation  field  differs  significantly  because  of  their  relative 
temperatures  and  proximity  to  the  nearly  black  underlying  surfaces.  Brightness  temperatures  above  the 
ICP  layer  converge  more  rapidly  with  increasing  optical  depth  than  do  corresponding  temperatures  above 
cirrus  because  the  top  of  the  ICP  layer  coincides  with  the  warmest  region  of  the  atmosphere  which  is 
directly  influenced  by  surface  emissions.  Radiation  emitted  by  the  surface  contributes  significantly  to  the 
total  upwelling  radiation  through  absorption  and  secondary  emission  by  the  ICP  layer  at  a  relatively  warm 
layer  temperature.  A  similar  radiative  effect  occurs  in  the  case  of  cirrus,  but  the  cold,  dry  atmosphere 
below  the  cirrus  layer  has  little  effect  on  the  upward  radiative  flux  and  the  cloud  particles  themselves 
absorb  and  re-emit  this  radiation  at  a  much  colder  temperature.  As  with  haze  and  stratiform  water  clouds, 
these  different  radiative  effects  give  rise  to  distinct  signatures  of  brightness  temperatures  and  bispectral 
differences  as  a  function  of  optical  depth. 
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II.2.4.2  NIR-IR  Bispectral  Differences 

Differences  in  BTD  signatures  for  combinations  of  surface  and  layer  types  as  a  function  of  optical 
depth  may  be  used  to  distinguish  varying  surface  and  atmospheric  properties  within  a  satellite  scene,  a 
necessary  step  in  developing  a  lead  detection  algorithm.  For  example,  at  a  cirrus  optical  depth  of  about 
3.0  the  BTD  between  channel  3  and  channel  4  (Tb3  -  T^)  over  an  open  lead  is  approximately  13K  whereas 
this  difference  is  only  about  2K  if  a  low  level  ICP  layer  of  equal  optical  depth  is  present  as  viewed  at 
nadir. 

To  more  clearly  illustrate  the  potential  use  of  split  window  imagery  to  distinguish  between  various 
cloud  and  aerosol  layers  that  are  common  in  the  Arctic  Figure  II. 5  was  constructed.  Shown  are  values 
of  (Tb3  -  Tm)  for  three  different  layers  already  discussed,  dry  (70%  RH)  and  wet  (99.9%  RH)  haze  layers 
within  the  boundary  layer  and  ICP  in  the  lowest  kilometer  of  the  atmosphere,  and  in  addition,  results  are 
given  for  the  stratus  layer  described  above  and  shown  in  Figure  n.3.  Each  panel  includes  plots  of  the 
simulated  Tb3  -  TM  values  for  the  four  surface  types  assuming  a  scan  angle  of  0°.  It  is  clear  that  bispectral 
differences  are  sensitive  not  only  to  optical  depth  but  also  to  relative  humidity  and  the  phase  of  the 
particles.  Bispectral  signatures  for  dry  aerosols  are  insensitive  to  increasing  optical  depth  except  for  those 
related  to  multiyear  ice.  Under  saturated  hazy  conditions  there  is  a  distinct  monotonic  decline  in  BTDs 
with  increasing  x  over  all  surface  types  with  negative  values  observed  except  over  multiyear  ice  when  x 
<.  5.  Such  separations  between  multiyear  ice  signatures  and  other  surface  types  should  permit  better 
identification  of  background  pixels  needed  to  normalize  thermal  contrasts  for  the  purpose  of  distinguishing 
leads.  If  a  stratus  layer  is  present  within  the  warm  region  of  the  atmosphere  a  sharp  fall-off  to 
significantly  large  negative  values  occurs  in  the  range  of  optical  deptns  between  0  end  about  2.5  with  a 
converging  upward  signal  as  x  increases  further.  In  the  case  of  a  low-lying  ICP  layer  positive  differences 
of  0.5K  to  3.5K  peaking  between  optical  depths  of  0.5  and  1.5,  are  apparent  whereas  tor  a  similar  f cirrus) 
layer  placed  high  in  the  atmosphere  BTDs  as  large  as  13K  were  noted  for  -  »  2.5  (sec  Figure  £1.4).  In 
theory,  if  TOA  NIR  and  ER  radiances  can  be  measured  accurately  much  information  can  be  extracted  by 
analyzing  bispectral  images.  Unfortunately,  the  AVHRR  channel  3  data  is  reported  to  be  too  noisy  to  be 
useful  for  cloud  detection  at  cold  temperatures  ( Yamanouchi  et  al.,  1987),  but  hopefully,  future  spacebome 
radiometers  will  provide  data  of  sufficient  quality  to  resolve  the  signatures  described  here. 
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Fig.  U.5.  AVHRR  channel  3  minus  channel  4  brightness  temperatures  o.s  a  function  of  0.55  urn  optical  depth  for 
four  intervening  layer  types:  aerosol  layers  having  mean  relative  humidities  of  70%  and  99.9%,  respectively,  ice 
crystal  precipitation  and  stratus  cloud.  Brightness  temperature  differences  are  shown  for  the  four  surface  types 
described  in  the  text. 
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II.2.4.3  Normalized  Atmospheric  Contrast 

Normalized  atmospheric  contrast  is  a  wavelength  dependent  quantity,  varying  with  atmospheric 
opacity,  expressed  in  terms  of  the  brightness  temperature  difference  between  any  given  pixel  and  the 
background  scene  normalized  by  the  brightness  temperature  of  the  background  (in  this  case  taken  to  be 
the  2  m  thick  multiyear  ice): 


C(x,X)  = 


Tt(x,X)-Tb(x,X) 

Tb^X) 


(1) 


where  Tr  is  the  brightness  temperature  of  the  target  (lead)  and  rB  is  the  back  .ound  brightness 
temperature.  Hereafter,  in  all  expressions  involving  contrast  the  wavelength  dependence  is  assumed  and 
the  X  is  omitted.  Figure  II.7  shows  the  behavior  of  this  quantity  for  IR  channel  4  derived  from  the  results 
presented  in  Figure  13.4.  This  measure  of  contrast  provides  a  potentially  powerful  means  to  detect  leads 
using  thermal  imagery  both  for  daytime  and  nighttime  conditions.  Although,  in  most  cases.  Figure  II.6 
indicates  a  rapid  decrease  in  normalized  contrasts  as  optical  depth  increases,  realistic  cirrus,  ICP  or  haze 
optical  depths  are  generally  within  a  range  that  should  permit  the  resolution  of  leads  using  thermal 
contrasts  provided  that  radiances  can  be  measured  accurately,  sensor  field-of-view  is  small  relative  to  lead 
widths  and  full  use  of  ancillary  data  is  made.  Curry  et  al.  (1990)  for  instance,  measured  ICP  visible 
optical  depths  ranging  from  about  0.03  to  20,  but  in  five  out  of  seven  cases  x  was  less  than  about  5,  within 
a  range  in  which  thermal  contrast  should  be  measurable.  With  regard  to  cirrus,  even  for  large  extinction 
coefficients  their  optical  depths  are  limited  because  they  are  confined  to  regions  of  the  upper  troposphere 
bounded  above  by  the  tropopause.  As  an  example,  a  cirrus  cloud  having  a  large  visible  volume  extinction 
coefficients,  say  k„,  =  0.2,  7  km  thick  (AZ  =  7.0)  would  have  an  optical  depth  of  only  1.4,  where  x  = 
k„,AZ.  Such  large  extinctions  would  exist  only  for  cirrus  composed  of  very  large  ice  crystals  with 
proportionally  large  total  ice  .vater  contents  (IWCs)  (e.g..  Stone  et  al.,  1990)  which  occur  rarely  in 
extremely  cold  environments  (e.g..  Stone,  1992;  Platt  and  Harshvardhan.  1988;  Heymsfield  and  Platt, 
1984).  In  cases  involving  haze,  the  aerosol  loading  would  need  to  be  extreme  before  thermal  contrasts 
diminished  significantly.  There  is  no  observational  evidence  that  aerosol  layers  can  attain  optical  depths 
that  would  preclude  lead  detection  using  thermal  imagery.  For  instance,  mean  0.5  pm  optical  depths  for 
haze  layers  observed  over  Barrow,  Alaska,  even  when  nearly  saturated,  were  only  about  0.2  ( Mendonca 
et  al.,  1981)  and  the  maximum  0.5  pm  optical  depth  measured  during  what  has  been  described  as  a 
"megahaze"  event  was  about  0.7  ( Dutton  et  al.,  1989).  Of  course,  mixed  phase  haze  layers  containing 
large  concentrations  of  ice  crystals  will  tend  to  attenuate  thermal  radiation  much  like  pure  ICP  layers  do. 
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Fig.  n.6.  Differences  between  AVHRR  channels  3  and 
4  brightness  temperatures  as  a  function  of  boundary 
layer  haze  optical  depth  at  0.55pm  during  the  day 
(thick  lines)  and  at  night  (thin  lines)  over  each  of  the 
four  surface  types  described  in  the  text.  Satellite  scan 
angle  is  50^,  solar  zenith  angle  is  75°  and  solar  azimuth 
angle  is  30°. 


Finally,  channel  contrasts  for  the  380  m  thick  stratus  cloud  (computed  but  not  shown)  were  also 
analyzed  to  evaluate  whether  or  not  leads  are  detectable  during  winter  when  stratus  layers  tend  to  be 
optically  thin.  IR  channel  brightness  temperatures  for  the  four  surface  types  were  observed  to  converge 
in  a  similar  manner  as  was  noted  for  the  boundary  layer  ICP  (Figure  II.4),  thus  simulated  IR  contrasts 
under  the  influence  of  stratus  clouds  are  nearly  equal  to  those  computed  for  the  ice  crystal  precipitation 
Lyer  shown  in  Figure  II.7.  The  convergence  of  channel  3  brightness  temperatures  with  optical  depth  was 
found  to  be  even  more  rapid  than  for  the  IR  channels  so  that  NIR  contrast  values  diminish  more  quickly 
with  optical  depth.  Both  IR  and  NIR  contrasts  are  likely  to  be  below  measurable  threshold  values 
(discussed  below)  to  make  lead  detection  possible  when  stratus  layers  are  present  because  their  optical 
depths  are  typically  in  the  range  of  2  or  greater.  Our  simulations  indicate  that  normalized  contrasts  will 
be  very  difficult  to  resolve  at  this  level  of  opacity. 
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Fig.  n.7.  Normalized  AVHRR  channel  4  contrasts  between  three  lead  ice  thicknesses  (open  water,  5  cm  and  15  cm 
ice)  and  the  background  2  m  thick  ice  for  boundary  layer  haze,  ice  crystal  precipitation  and  high-level  cirrus  cloud 
as  a  function  of  0.55pm  optical  depth.  Results  are  shown  for  satellite  view  angles  of  0“  and  SIF  derived  from  the 
brightness  temperatures  shown  in  Figure  n.4. 
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II.2.4.4  Lead  Detection  Based  on  Critical  Contrast 

The  normalized  contrast  C(xX)  in  (1)  is  defined  in  terms  of  the  target  (lead)  and  background 
(multiyear  ice)  temperatures  but  says  nothing  about  the  geometrical  characteristics  of  the  target  or  imaging 
system.  In  such  a  context  it  assumes  that  a  given  lead  pixel  is  completely  within  the  FOV  of  the  satellite 
radiometer.  We  now  examine  how  contrast  varies  as  a  function  of  lead  width  and  sensor  field-of-view 
to  account  for  image  pixels  that  contain  both  lead  and  multiyear  ice  surface  types. 

Letting  p  be  the  fractional  area  covered  by  a  lead  within  a  pixel;  i.e.,  p  =  width/FOV,  the  total 
contrast  C,0,  that  takes  into  account  the  reduction  in  temperature  contrast  due  to  atmospheric  and  spatial 
effects  is 


Ctot(  T)  - 


[prr(x)+(i-p)re(x)]-re(T) 


=  pC(x) 


W 


If  every  pixel  in  the  image  is  to  be  labeled  as  either  a  lead  pixel  or  a  background  pixel,  then  some 
thresholding  operation  must  be  used.  One  possible  method  is  to  choose  as  a  threshold  the  mean 
background  temperature  plus  some  multiple  of  its  variability  a,  say  TB(x)+2o.  (In  reality  a  may  also  be 
a  function  of  x.)  As  in  (1),  this  threshold  can  also  be  expressed  as  a  normalized  (non-dimensional) 
contrast  ratio: 


If  the  observed  total  contrast  of  a  pixel  is  less  than  this  value,  then  the  pixel  is  not  a  lead  pixel.  This 
threshold  contrast  includes  implicitly  the  effect  of  the  fractional  area  covered  by  a  lead  within  a  pixel;  i.e., 
it  is  a  total  contrast.  Low  values  of  y  will  generally  result  in  more  pixels  being  labelled  as  lead  pixels 
because  the  background  is  more  homogeneous  if  o  is  small.  Thermal  features  should  be  more 
distinguishable  when  contrasted  with  a  more  homogeneous  background  scene,  y  can  also  be  defined  in 
terms  of  some  critical  normalized  atmospheric  contrast: 

C*(x,y)  =  1 

P 

where  the  asterisk  represents  a  critical  (cutoff)  value.  This  expresses  the  normalized  atmospheric  contrast 
necessary  for  a  lead  to  be  detected  if  an  intervening  layer  of  optical  depth  x  is  present. 

To  address  the  question  of  what  minimum  lead  width  can  actually  be  resolved  under  specified 
atmospheric  conditions  and  sensor  FOV,  we  need  to  eliminate  hypothetically  the  atmospheric  effects  and 
account  for  varying  FOV.  We  therefore  define  the  critical  contrast  of  a  lead  as  the  normalized 
atmospheric  contrast  at  zero  optical  depth  C*(x=0,y).  We  can  relate  C‘(x,y)  back  to  the  critical  contrast 
C*(x=0,y)  with  the  data  provided  in  Figure  II.7  in  the  following  manner.  Nadir  viewing  is  assumed.  A 
log  surface  is  fitted  to  the  contrast  data  for  ICP  and  cirrus  cloud,  while  a  planar  surface  is  fitted  to  the 
aerosol  contrast  data.  Least  squares  regression  is  employed  where  the  independent  variables  are  the  actual 
brightness  temperatures  at  x=0  of  the  three  lead  types  and  optical  depth.  The  dependent  variable  is 
contrast: 
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where  the  subscript  i  denotes  aerosol,  ICP  or  cirrus  cloud;  C*'i=ln(C')  for  ICP  or  cirrus  cloud  and  C'~C 
for  an  aerosol  layer.  For  a  given  width  to  FOV  ratio  p,  optical  depth  x  and  threshold  contrast  y,  the  target 
brightness  temperature  Tr  is  estimated  for  the  appropriate  background  temperature  from  Figure  n.4.  Tr 
is  then  used  in  the  regression  equation  above  to  determine  the  atmospheric  contrast  at  that  optical  depth. 
Inverting  the  equation  and  setting  x  equal  to  0  gives  the  critical  contrast  C* (T=0,y).  As  discussed  earlier, 
IR  brightness  temperatures  at  t=0  are  very  close  to  the  physical  surface  temperatures  because  these 
surfaces  have  emissivities  approaching  unity  and  the  low  water  vapor  content  of  the  atmosphere  has  a 
negligible  influence  on  the  upwelling  radiation. 

Figure  H.8  shows  critical  normalized  contrast  contours  for  leads  as  a  function  of  p  and  optical 
depth  for  haze,  ICP  and  cirrus  doud  at  two  assumed  values  of  y.  As  mentioned  above,  it  is  unlikely  that 
contrasts  in  the  presence  of  stratus  clouds  can  be  sufficiently  resolved  for  realistic  values  of  stratus  optical 
depth  to  permit  lead  detection,  therefore  stratus  is  omitted  from  this  analysis.  The  contours  in  Figure  II.8 
indicate  the  critical  contrast  that  a  pixel  must  exceed  in  order  to  be  designated  as  a  lead  pixel  for  the 
assumed  threshold  contrast  y,  again  assuming  mean  clear-sky  January  conditions  in  the  central  Arctic. 
Such  plots  can  be  used  to  estimate  the  minimum  lead  width  resolvable  in  an  AVHRR  channel  4  image 
under  certain  conditions  in  the  following  manner.  Assume,  for  example,  that  the  sensor’s  resolution  is 
1.0  km  (FOV  =  1.0)  at  nadir  and  a  1  km  thick  ICP  layer  is  present  above  the  surface.  If  we  prescribe  a 
normalized  contrast  threshold  of  0.04  and  a  critical  contrast  of  0.2  as  detection  criteria,  then  Figure  II.8 
can  be  used  to  estimate  the  narrowest  resolvable  lead  in  a  channel  4  AVHRR  image  as  a  function  of  the 
layer’s  optical  depth.  For  an  optical  depth  of  8,  the  width/FOV  ratio  is  0.35  thus,  the  narrowest  detectable 
lead  would  be  0.35  km  wide.  Under  these  conditions,  any  L  id  less  than  about  350  m  wide  would  not  be 
detectable  if  the  ICP  optical  depth  exceeded  8.0.  Following  the  same  procedure  but  assuming  an  ICP 
optical  depth  of  1,  the  minimum  detectable  lead  width  is  about  0.22  km  given  equal  threshold  criteria. 
If  we  now  relax  the  threshold  criteria  and  do  the  same  analysis  for  y  =  0.02,  then  a  lead  would  need  to 
be  only  about  180  m  wide  in  order  to  be  resolved  at  a  pixel  resolution  of  1  km  if  x  =  8,  or  1 10  m  wide 
if  t  =  1.  For  any  given  combination  of  optical  depth,  contrast  threshold  and  critical  contrast  the  minimum 
detectable  lead  width  will  be  systematically  smaller  if,  instead  of  ICP,  a  haze  layer  is  present  and  slightly 
greater  if  cirrus  cloud  is  present  based  on  this  technique  and  the  theoretical  results  presented  in  Figure  DL8. 

It  is  suggested  here  that  bispectral  techniques  be  used  in  conjunction  with  image  contrast  analyses 
to  develop  an  operational  procedure  to  detect  and  map  leads  in  the  polar  regions.  However,  because  a 
continuum  of  signatures  exist  depending  on  atmospheric,  surface  and  geometric  effects,  it  will  be  essential 
to  constrain  the  problem  by  first  determining  the  state  of  the  atmosphere  and  surface  using  a  combination 
of  multispectral  techniques.  Of  particular  importance  will  be  the  retrieval  of  surface  temperatures  (e.g.. 
Key  and  Haefliger,  1992)  because  the  thermal  contrast  between  leads  and  the  background  ice  pack  may 
be  a  key  parameter  for  determining  the  threshold  of  lead  detection  under  various  conditions. 
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Fig.  E8.  Critical  contrast  of  a  lead  as  a  function  of  its  fractional  coverage  within  an  image  pixel  (width/FOV)  and 
the  0.55pm  optical  depth  of  aerosols  (haze),  ice  crystal  precipitation  and  cirrus  cloud.  For  a  given  set  of  geometrical 
and  atmospheric  conditions  the  contours  indicate  the  minimum  contrast  (at  an  optical  depth  of  zero)  that  must  exist 
in  order  for  a  lead  to  be  detectable  for  specified  values  of  threshold  contr_  t  y  assuming  a  satellite  viewing  angle  of 
0°.  Mean  clear-sky  January  conditions  for  the  central  Arctic  are  assumed.  Results  are  given  for  two  values  of 
threshold  contrast:  y=0.04  and  y=0.02  as  defined  in  the  text. 
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II.2.5  DISCUSSION 

The  theoretical  results  presented  above  are  enlightening,  but  we  have  obviously  oversimplified  the 
problems  associated  with  lead  detection.  We  have  presented  results  that  only  represent  average  January 
conditions  in  the  central  Arctic  and  consider  only  four  discrete  surface  types  and  four  idealized 
hypothetical  models  to  simulate  what  in  nature  is  a  complicated  mix  of  intervening  atmospheric  effects. 
In  reality,  any  of  these  variables  assume  a  continuum  of  values  that  change  spatially  and  temporally. 
Clouds  and  aerosol  layers  are  naturally  inhomogeneous  having  physical,  radiative  and  microphysical 
properties  that  vary  in  space  and  time  and  frequently  occur  as  multiple  layers  of  mixed  phase  particles. 
Intense  winds,  for  example,  dynamically  force  stratiform  layers  sometimes  creating  banded  structures, 
especially  for  layers  composed  of  condensed  particles  downwind  of  leads.  The  detection  of  leads  is 
further  complicated  by  constantly  changing  viewing  geometries  related  to  sensor  FOV,  satellite  scan  angle 
(and  sun  angle  if  daytime  conditions  exist).  Non-linear  radiative  effects  are  caused  by  increasing  optical 
paths  at  the  same  time  that  pixel  resolution  degrades  with  increasing  scan  angle.  To  develop  an 
operational  lead  detection  algorithm,  highly  parameterized  models  will  need  to  be  used,  perhaps  in 
conjunction  with  comprehensive  "look-up"  tables  listing  expected  contrast  values  and  corresponding  BTDs 
for  a  realistic  range  of  combined  surface  and  atmospheric  properties.  Furthermore,  a  step-wise  approach 
will  be  necessary  utilizing  ancillary  information  to  further  constrain  the  problem.  The  use  of  cloud 
clearing  algorithms  is  essential  to  assure  accurate  surface  temperature  retrievals  and  remote  sounding 
techniques  need  to  be  improved  in  order  to  resolve  atmospheric  temperature  and  humidity  profiles  and  to 
determine  the  physical,  radiative  and  microphysical  properties  of  intervening  layers.  Theoretically, 
multispectral  analyses  provide  a  basis  for  estimating  these  properties,  but  in  reality  it  may  be  years  before 
prototype  methods  are  validL.ed  and  approved  for  operational  use. 
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The  ability  to  retrieve  surface  parameters  from  satellite  data  in  the  polar  regions  is  constrained  by 
our  limited  knowledge  of  atmospheric  temperature,  humidity,  and  aerosol  profiles,  the  microphysical 
properties  of  polar  clouds,  and  the  spectral  characteristics  of  the  wide  variety  of  surface  types  found  there. 
In  this  section  we  present  results  in  the  retrieval  of  ice  surface  temperature  (1ST)  from  the  thermal 
channels  of  the  Advanced  Very  High  Resolution  Radiometer  (AVHRR)  sensors  on-board  the  NOAA  series 
satellites. 

Sea  and  land  surface  temperature  (SST  and  LST)  retrieval  algorithms  have  been  developed  using 
the  thermal  infrared  window  portion  of  the  spectrum,  with  the  degree  of  success  dependent  primarily  upon 
the  variability  of  the  surface  and  atmospheric  characteristics.  The  general  approach  to  estimating  surface 
temperature  is  to  relate  satellite  observations  to  surface  temperature  observations  with  a  regression  model. 
Lacking  sufficient  observations,  however,  satellite  radiances  or  brightness  temperatures  can  be  modeled 
by  application  of  the  radiative  transfer  equation.  This  approach  is  commonly  used  for  SST  retrieval. 

To  our  knowledge,  little  effort  has  been  directed  to  the  retrieval  of  the  sea  ice  surface  temperature 
(1ST)  in  the  arctic,  an  area  where  the  first  effects  of  a  changing  climate  are  expected  to  be  seen.  The 
reason  is  not  one  of  methodology,  but  rather  our  limited  knowledge  of  atmospheric  temperature,  humidity, 
and  aerosol  profiles,  the  microphysical  properties  of  polar  clouds,  and  the  spectral  characteristics  of  the 
wide  variety  of  surface  types  found  there.  We  have  developed  a  means  to  correct  for  the  atmospheric 
attenuation  of  satellite-measured  clear  sky  brightness  temperatures  used  in  the  retrieval  of  ice  surface 
temperature  from  the  split-window  thermal  channels  of  the  AVHRR  sensors  on-board  three  of  the  NOAA 
series  satellites.  These  corrections  are  specified  for  three  different  "seasons"  and  as  a  function  of  satellite 
viewing  angle,  and  are  expected  to  be  applicable  to  the  perennial  ice  pack  in  the  central  Arctic  Basin 
(Figure  II.  1).  We  do  not  develop  a  completely  new  methodology;  instead  we  modify  a  standard  procedure 
for  use  with  arctic-specific  data.  It  is  assumed  that  a  valid  cloud-clearing  algorithm  exists  and  that  only 
clear  sky  radiances  are  being  examined.  The  cloud  clearing  problem  in  polar  satellite  data  is  not  trivial, 
however.  For  a  review  of  polar  cloud  detection  algorithms,  see  Key  and  Barry  (1989)  and  Sakellariou 
et  al.  (1991). 

For  the  retrieval  of  1ST  a  multi-channel  algorithm  that  uses  empirical  relationships  to  correct  for 
water  vapor  absorption  is  employed; 

Tice  -  a  +  bT 4  +  cT5  +  d[(r4-r5)sec0] 

where  T4  and  Ts  are  the  satellite-measured  brightness  temperatures  (K)  in  the  AVHRR  thermal  channels 
and  0  is  the  sensor  scan  angle.  The  coefficients  are  determined  through  a  least  squares  regression 
procedure,  where  surface  temperatures  are  regressed  against  modeled  brightness  temperatures. 

AVHRR  thermal  channel  radiances  are  simulated  with  LOWTRAN  as  described  previously. 
Directional  surface  emissivities  for  snow  are  modeled  ( Dozier  and  Warren,  1982):  the  single  scattering 
albedo  and  asymmetry  factor  in  the  scattering  phase  function  are  calculated  from  the  Mie  equations  and 
the  directional,  wavelength-dependent  emissivities  are  derived  from  the  delta-Eddington  approximation  to 
the  equation  of  radiative  transfer.  The  directional  emissivities  are  then  integrated  with  the  response 
function  for  each  channel  as  described  previously. 

The  use  of  the  rawinsonde  profiles  in  modeling  the  surface  temr,,’rature  requires  an  additional  step 
since  the  first  measurement  in  each  profile  is  the  shelter  temper  .  not  the  surface  temperature. 
Therefore,  the  (unknown)  surface  temperature  for  each  profile  is  aligned  a  series  of  values  representing 
the  range  of  possible  surface  temperatures  for  the  observed  conditions  during  the  month  to  which  the 
rofile  belongs.  An  energy  balance  model  ( Maykut ,  1982)  is  used  to  determine  these  surface  temperatures, 
based  on  the  observed  range  of  shelter  temperatures  and  wind  speeds  (the  mean  ±1  standard  deviation) 
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in  the  ice  island  data  for  each  month. 

The  seasonal  dependence  of  the  coefficients  is  illustrated  in  Table  II.4,  where  coefficients  from 
each  season  were  applied  to  data  from  every  other  season.  Results  are  shown  for  NOAA  9  and  indicate 
errors  between  0. 1  K  for  transition  coefficients  with  winter  data  and  0.6  K  when  summer  coefficients  are 
used  with  winter  data.  Similarly,  the  satellite  dependence  of  the  coefficients  is  shown  in  Table  II.5  for 
summer  conditions.  On  the  average,  errors  ranging  from  0.1  to  1.0  K,  depending  on  season,  can  be 
expected  when  applying  coefficients  derived  for  one  satellite  to  data  from  another,  the  smallest  errors 
occurring  between  NOAA  7  and  9  coefficients  and  data.  Using  SST  coefficients  developed  for  the  North 
Atlantic  and  the  Greenland  Sea  area  to  estimate  1ST  would  result  in  an  underestimate  of  up  to  0.7  K, 
largest  in  winter  and  at  scan  angles  of  40°  and  greater.  While  the  sensor  scan  angle  is  included  explicitly 
in  the  correction  equation,  its  effect  in  the  dry  arctic  atmosphere  is  small,  generally  less  than  0. 1  K. 

Surface  temperature  measurements  taken  by  a  PRT-5  thermal  radiometer  during  CEAREX  in 
March  1989  were  compared  to  estimated  ISTs  from  NOAA  1 1  AVHRR  data.  The  mean  1ST  for  a  sample 
of  four  AVHRR  pixels  was  258.9  K  while  the  mean  PRT-5  temperature  (adjusted  for  an  emissivity  of 
0.998)  of  four  consecutive  measurements  one  kilometer  apart  was  259.04  K.  Given  the  difficulties  in 
comparing  the  two  data  sets  these  results  are  encouraging. 

In  summary,  using  the  split  window  channels  and  scan  angle,  the  rms  error  in  the  estimated  ice 
surface  temperature  is  less  than  0.1  K  in  all  seasons.  Inclusion  of  channel  3  (3.7  pm)  during  the  winter 
decreases  the  rms  error  by  less  than  0.003  K.  Overall,  employing  the  1ST  coefficients  results  in  increased 
accuracy  of  up  to  0.6  K  over  SST  coefficients  developed  for  the  North  Atlantic  and  the  Greenland  Sea 
areas. 

Additional  details  are  provided  in  Key  and  Haefliger  (1992). 


Table  II.4.  RMS  error  in  applying  coefficients  (NOAA  9)  developed  for  one  season  (left)  to  data  from 
another  (top). 


Coefficients: 

Data  from: 
Winter  Summer 

Transition 

Winter 

0 

0.403 

0.128 

Summer 

0.587 

0 

0.342 

Transition 

0.117 

0.219 

0 
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Table  H.5.  RMS  error  applying  coefficients  (summer)  developed  for  one  satellite  (left) 
to  data  from  another  (top). 


Data  from: 

NOAA7 

NOAA  9 

NOAA  11 

Coefficients: 

NOAA 7  0 

0.272 

0.655 

NOAA  9  0.296 

0 

1.017 

NOAA  11  0.682 

0.96! 

0 
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Part  iii:  Analysis  of  Leadex  and  Other  In  Situ  Data 


In  this  section  we  discuss  data  collected  during  LEADEX  and  its  application  to  some  of  the 
procedures  and  models  described  in  the  previous  sections.  These  applications  are  not  exhaustive,  but  do 
provide  some  insight  into  problems  that  remain  in  the  satellite  remote  sensing  of  leads. 

m.l  ARCTIC  AEROSOLS 

The  June  1991  eruptions  of  Mount  Pinatubo  ejected  massive  amounts  of  debris  and  sulfur  dioxide 
gas  into  the  stratosphere.  Increased  concentrations  of  aerosols  in  the  stratosphere  can  perturb  the  radiative 
balance  of  the  entire  earth/atmosphere  system  (e.g.,  Minnis  et  al.,  1993;  Hansen  et  al.,  1992;  Dutton  and 
Christy,  1992).  This  "volcanic  forcing"  is  characterized  by  tropospheric  cooling  caused  by  an  increase 
in  the  planetary  albedo  and  by  stratospheric  wanning  primarily  due  to  infrared  absorption  by  aerosol 
particles.  After  reach'  'g  a  peak,  the  stratosphere’s  opacity  normally  decays  exponentially  in  time  at  a  rate 
dependent  on  the  magnitude,  time  of  year  and  location  of  the  eruption  (e.g.,  Gerber  and  Deepak,  1984; 
Hofmann  and  Rose::,  1987).  Although  the  net  radiative  effect  of  atmospheric  aerosols  depends  on  many 
factors,  the  optical  thickness  and  effective  size  distribution  of  the  aerosols  are  most  important  ( Lads  et 
al.,  1992).  Aerosol  optical  thickness  is  a  nondimensional  parameter  used  to  quantify  the  spectral 
extinction  of  direct  solar  irradiance  by  aerosols  integrated  along  a  path  between  an  observer  and  the  sun. 
If  measurements  of  extinction  can  be  obtained  over  a  suitable  range  of  wavelengths,  then  an  effective  size 
distribution  of  an  aerosol  layer  can  be  inferred  (King  et  al.,  1978). 

During  the  spring  of  1992  an  extensive  series  of  in  situ  measurements  were  made  using  airborne 
techniques  as  pan  of  the  Fourth  Arctic  Gas  and  Aerosol  Sampling  Program  (AGASP-IV)  in  conjunction 
with  the  Arctic  Leads  Experiment  (LEADEX).  Nearly  1300  spectral  measurements  of  solar  irradiance 
were  made  from  near  the  surface  into  the  stratosphere  using  handheld  sunphotometers  during  seven  flights 
of  the  NOAA  WP-3D  aircraft.  We  focus  here  on  an  analysis  of  the  stratospheric  data  to  quantify  the 
spectral  opacity  and  infer  effective  size  distributions  for  the  Pinatubo  aerosols  present  in  the  Arctic. 
Ancillary  surface  measurements  are  presented  in  support  of  the  aircraft  data  analysis  and  are  further  used 
to  estimate  a  decay  rate  of  stratospheric  optical  depth  following  the  period  of  peak  aerosol  concentration. 

While  the  radiative  transfer  modeling  results  presented  in  previous  sections  indicates  that  aerosols 
play  a  relatively  small  role  in  the  remote  sensing  of  leads  using  thermal  data,  the  results  presented  in  this 
section  are  important  because 

1.  the  aerosol  amounts  measured  are  twice  the  value  typical  of  the  Arctic  so  that  their  impact  on 
leads  remote  sensing  in  the  thermal  portion  of  the  spectrum  cannot  be  ignored,  and 

2.  multiple  scattering  of  shortwave  radiation  by  aerosols  is  significant,  especially  at  the  levels 
measured. 

While  we  have  not  emphasized  remote  sensing  with  shortwave  sensor  channels,  they  provide  a  potentially 
valuable  source  of  information  and  should  be  the  focus  of  future  research. 
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III.  1.1  MEASUREMENTS 

The  data  are  derived  from  two  types  of  radiometric  observations:  (1)  airborne  and  surface-based 
measurements  made  using  sunphotometer  s,  and  (2)  direct  beam,  wideband  solar  irradiance  measurements 
made  at  a  ground  station.  The  primary  data  are  derived  from  spectral  measurements  made  during  flight 
segments  flown  above  the  tropopause.  Additional  sunphotometer  measurements  were  made  at  the  surface 
near  Resolute,  N.W.T.,  as  part  of  the  1992  and  1993  Seasonal  Sea  Ice  Monitoring  and  Modeling  Site 
(SIMMS)  field  programs  ( Reddan  et  al.,  1992),  and  in  the  vicinity  of  Anchorage,  Alaska.  The 
sunphotometer  observations  were  made  using  two  handheld,  dual-channel  instruments  that  sense  directly 
transmitted  solar  irradiance  at  380  and  500  nm,  and  at  778  and  862  nm,  respectively;  each  channel  having 
a  nominal  half-bandwidth  of  5  nm  and  a  field  of  view  of  2.4°.  The  wideband  (350-695  nm)  pyrheliometer 
data  were  collected  at  the  NOAA  Climate  Monitoring  and  Diagnostics  Laboratory’s  Barrow  Observatory 
(CMDL/BRW).  The  "wideband  method"  used  to  estimate  optical  depth  was  described  by  Dutton  and 
Christy  (1992).  Only  data  collected  during  cloud-free  periods  are  analyzed.  The  locations  and  dates 
corresponding  to  the  various  measurement  periods  are  shown  in  Figure  III.  1.  The  curved  vectors  are  back- 
trajectories  representing  stratospheric  winds  prior  to  each  flight. 


Fig.  m.l.  Distribution  of  sunphotometer  measurements  made  during  AGASP-IV/LEADEX  stratospheric  flight 
segments  and  at  surface  locations  (ANCH,  TALK,  and  SIMMS).  Wideband  pyrheliometric  measurements  were  made 
at  the  Barrow  Observatory  (BRW).  Curved  vectors  represent  stratospheric  winds  for  36  hours  prior  to  each  flight 
(back-trajectories). 
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Airborne  observations  were  made  through  optical  glass  windows  whenever  viewing  was  possible 
within  ±30°  of  the  solar  azimuth  relative  to  a  plane  perpendicular  to  the  aircraft  heading  (e.g.,  Dutton  et 
al.,  1989).  The  sunphotometers  were  carefully  calibrated  at  high-altitude  sites  before  °nd  after  AGASP-IV 
and  again  prior  to  SIMMS  ’93  in  accordance  with  recommended  procedures  using  the  Langley  plot  method 
(Shaw,  1983).  Each  instrument  was  determined  to  be  stable  with  precision  of  about  0.002  optical  depth 
units  (ODUs).  All  optical  depth  values  reported  in  this  paper  account  for  Rayleigh  scattering,  ozone 
absorption  at  500  nm,  and  changes  in  relative  airmass  as  a  function  of  time  and  location.  Attenuation  by 
the  optical  glass  was  also  accounted  for  when  reducing  the  aircraft  data.  Thus,  from  total  slant-path 
irradiance  measurements,  columnar  aerosol  optical  depths  were  computed.  Estimated  uncertainties  due 
to  all  sources  of  error  are  no  more  than  0.02  ODUs  (Reddy  et  al,  1990),  or  less  than  10%  of  typical 
values  reported  here.  The  wideband  measurements  are  accurate  to  within  ±0.04  ODUs  (Dutton  and 
Christy,  1992). 

III.  1.2  RESULTS 

III.  1.2.1  Stratr  spheric  Aerosol  Optical  Depths 

Ln.  y  direct  measurements  made  above  the  tropopause  and  the  clear-sky,  ground-based 
measurements  were  selected  for  analysis.  Tropopause  heights  for  each  flight  were  determined  on  the  basis 
of  the  analyses  of  Herbert  et  al.  (1993);  these  range  from  about  7  km  to  10  km,  depending  on 
geographical  location  and  synoptic  conditions.  Stratospheric  components  were  estimated  from  the  surface 
measurements  of  total-column  optical  depth  by  subtracting  predetermined  "background"  values  for  the 
troposphere,  a  method  also  used  by  Asano  et  al.  (1993).  For  our  purposes  the  March  1979-1982  monthly 
mean  background  values  of  spectral  aerosol  optical  depth  estimated  for  BRW  (Dutton  et  al.,  1984)  were 
systematically  subtracted  from  the  respective  measured  values.  Note,  however,  that  March  background 
values  are  typically  smaller  than  those  for  April  and  larger  than  those  for  May  or  Tune  based  on  recent 
wideband  analyses  (Dutton,  unpublished  data).  Figure  III.2  shows  the  mean  stratospheric  aerosol  optical 
depths,  plus  and  minus  one  standard  deviation  (±lo),  derived  from  observations  made  at  the  corresponding 
times  and  locations  shown  in  Figure  III.  1 .  The  tropospheric  background  values  used  to  correct  the  surface 
measurements  are  also  shown  with  ranges  of  uncertainty  indicated. 

Several  features  of  Figure  III.2  are  notable.  First  and  most  important,  the  values  derived  for  all 
times  and  locations  are  1  to  2  orders  of  magnitude  greater  than  stratospheric  background  levels  (Toon  and 
Pollack,  1976);  these  values  also  exceed  similar  measurements  made  in  the  Arctic  approximately  a  year 
after  El  Chichon  erupted  (Spinhirne  and  King,  1985;  Dutton  et  al.,  1984).  Second,  optical  depths  derived 
from  the  aircraft  observations  tend  to  be  rather  flat  spectrally  in  the  visible  range  compared  with  the  1992 
surface  results.  Third,  three  of  the  flights,  403,  406  and  407,  indicate  a  high  degree  of  homogeneity  in 
time  and  space,  as  evidenced  by  their  small  standard  deviations,  similar  magnitudes,  and  spectral 
dependencies.  Fourth,  the  values  for  flight  402  are  about  60%  lower  than  those  measured  during  flight 
407  despite  their  having  similar  flight  tracks  and  altitudes  relative  to  the  tropopause.  Fifth,  the  curve  for 
SIMMS '92  also  falls  below  all  of  the  AGASP  curves  except  for  that  of  flight  402,  suggesting  that  either 
temporal  and/or  spatial  variations  occurred  over  the  period  and  geographical  region  represented  by  these 
data,  or  that  the  SIMMS  data  are  negatively  biased  due  to  incorrect  assumptions  made  regarding  the 
tropospheric  background  corrections.  And  last,  the  values  for  SIMMS '93  are  significantly  lower  and  have 
a  different  spectral  signature  than  those  for  SIMMS '92,  indicating  a  decay  in  opacity  from  one  year  to  the 
next  as  well  as  a  change  in  microphysical  characteristics. 
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Fig.  EL2.  Mean  stratospheric  spectral  aerosol  optical  depths  derived  for  the  flight  segments  and  surface  locations 
shown  in  Figure  DI.1.  Slight  spectral  offsets  are  used  to  clearly  indicate  values  at  plus  and  minus  one  standard 
deviation  (vertical  bars).  S92  and  S93  are  averages  of  the  SIMM’92  and  SIMMS’93  data,  respectively  (only  380 
and  500  nm  measurements  were  made  in  1993);  ANCH  is  the  average  of  the  Anchorage  (ANCH)  and  Talkeetna 
(TALK)  data.  The  TBG  curve  (adapted  from  Dutton  et  al.  (1984))  gives  the  values  and  estimated  ranges  of  the 
tropospheric  background  corrections  used  to  derive  stratospheric  optical  depths  from  the  surface-based  measurements. 


III.  1.2.2  Inferred  Aerosol  Size  Distributions 

Based  on  the  data  presented  in  Figure  M.2,  effective  aerosol  size  distributions  were  inferred  using 
the  constrained  linear  inversion  algorithm  of  King  et  al.  (1978).  The  radius  sensitivity  (rm,„  £  r  <.  r^ 
( Spinhirne  and  King,  1985)  determined  for  our  particular  set  of  measurements  was  within  the  range 
rm,n=0-  ±0.02  pm  and  ^=1.10  ±0.10  pm.  We  assumed  an  index  of  refraction  of  1.45-0.0/  based  on 

earlier  in  situ  observations  of  the  Pinatubo  aerosol  layers  ( Deshler  et  al.,  1992).  The  inversion  results  are 
presented  in  Figure  III.3.  Each  curve  shows  the  total  number  concentration  dN  (cm-2)  for  seven  radius 
increments  [dlog(r)].  The  vertical  bars  indicate  the  range  of  number  density  determined  by  inverting  the 
mean  spectral  optical  depth  data  ±la.  We  find  that,  for  the  period  and  geographical  region  of  interest, 
the  inferred  size  spectra  tend  to  fall  into  two  groups.  Both  are  bimodal,  having  a  large-particle  mode 
centered  at  about  0.50  pm  radius  and  a  small-particle  mode  of  higher  concentration  peaking  below  about 
0.18  pm.  These  results  suggest  that  the  volcanic  aerosols  present  in  the  Arctic  about  10  months  after 
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Pinatubo  erupted  were  slightly  larger  than  the  newly  formed  particles  observed  over  Laramie,  Wyoming 
several  weeks  after  the  eruption  ( Deshler  et  al.,  1992)  but  were  smaller  than  those  estimated  by  Asano  et 
al.  (1993)  in  Japan  approximately  2  months  before  the  Arctic  observations  were  made.  Their  independent 
determinations  of  size  spectra  also  show  the  bimodal  feature  noted  here  attributed  to  the  superposition  of 
a  monodisperse  large-particle  (volcanic)  mode  onto  a  small-particle  (background)  mode.  Similar  bimodal 
size  spectra  were  inferred  using  photometric  measurements  made  at  high  northern  latitudes  about  a  year 
after  El  Chichon  erupted  ( Spinhirne  and  King,  1985). 


Fig.  m.3.  Effective  aerosol  size  distributions  showing  the  range  of  number  concentrations  dN  (cm  21  at  each  radius 
interval  [dlog(r)]  inferred  from  the  optical  depth  data  shown  in  Figure  IH2. 


III.  1.2.3  Time  Decay  of  Stratospheric  Opacity 

To  evaluate  the  decay  of  the  Pinatubo  aerosol  layer(s),  we  analyzed  the  BRW  optical  depth  data 
and  the  spectrally-weighted  average  of  the  380  and  500  nm  SIMMS  observations  for  the  successive  1992 
and  1993  spring  periods.  At  BRW  the  optical  depth  reached  an  average  peak  value  of  about  0.19  during 
early  May  1992  and  declined  thereafter  (updated  from  Dutton  and  Christy  (1992)).  Assuming  that  optical 
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depth  t  decays  exponentially  in  time  t  after  reaching  its  peak  value  x  k,  i.e.,  i(t)  =  'tpeakexp(-//7'),  a  fit  to 
the  BRW  data  yields  a  time  constant  T  of  13.3  months.  A  corresponding  fit  to  the  SIMMS  data  yields 
a  time  constant  of  13.7  months,  assuming  that  the  peak  opacity  also  occurred  during  May  1992  and  that 
the  BRW  seasonal  background  tropospheric  corrections  were  valid  for  the  SIMMS  site.  We  give  credence 
to  these  estimates  because  independent  measurements  from  geographically  distinct  sites  were  used  to 
obtain  very  similar  results  noting,  however,  that  the  actual  overall  decay  of  volcanic  aerosols  may  not  be 
well  represented  by  a  smooth  exponential  function  ( Hofmann  and  Deshler,  1987).  At  least  another  year 
of  monitoring  is  needed  before  our  results  can  be  corroborated.  Our  preliminary  analysis  suggests  that 
the  Pinatubo  aerosols  are  being  systematically  removed  from  the  Arctic  stratosphere,  but  at  a  slower  rate 
than  estimated  following  earlier  volcanic  eruptions.  On  the  basis  of  in  situ  midlatitude  measurements  of 
aerosol  size  distributions,  Hofinann  and  Deshler  (1987)  estimated  the  decay  rate  (e-folding  time)  of  total 
stratospheric  mass  after  the  El  Chichon  eruption  to  be  about  10.3  months  and  between  8  and  10  months 
following  the  1974  eruption  of  Fuego  in  Guatemala.  High-latitude  satellite  observations  of  1.0  pm 
stratospheric  aerosol  optical  depth  ( McCormick  and  Trepte,  1987)  also  exhibited  faster  decay  rates 
following  the  1980  Mt.  St.  Helens  and  the  El  Chichon  post-volcanic  periods  than  we  report  here. 

III.  1.3  DISCUSSION 

It  appears  from  Figure  III.2  that  temporal  and/or  spatial  variations  in  stratospheric  aerosols 
occurred  in  the  Arctic  during  spring  1992.  The  lower  relative  opacity  noted  for  flight  402  can  be 
explained  by  differences  in  synoptic  conditions  during  the  respective  flight  periods  ( Herbert  et  al. ,  1993). 
During  flight  402  strong  northerly  winds  transported  polar  air  into  the  region  whereas  for  the  later  three 
flights  moderately  weak  southerly  flow  was  generally  observed.  An  analysis  of  isentropic  back-trajectories 
(e.g.,  Harris  and  Bodhaine,  1983)  based  on  the  ECMWF  2.5°  gridded  data  supports  the  synoptic  analyses. 
Figure  III.  1  shows  36  hour  back-trajectories  representing  the  flow  at  altitudes  ranging  from  about  14.5  to 
15.8  km  referenced  to  the  respective  midflight  segments.  This  analysis  suggests  that  relatively  clean 
Arctic  air  displaced  or  mixed  with  lower  latitude  air,  effectively  reducing  the  stratosphere’s  opacity  prior 
to  flight  402.  The  formation  of  the  polar  vortex  during  the  previous  autumn/winter  probably  prevented 
high  concentrations  of  Pinatubo  aerosols  from  penetrating  the  central  Arctic.  Similar  gradients  related  to 
the  position,  size,  and  shape  of  the  polar  vortex  were  observed  several  months  after  El  Chichon  erupted 
C McCormick  et  al.,  1983). 

The  relatively  large  380  nm  optical  depths  measured  during  flight  402  and  at  all  the  surface 
locations  during  spring  1992  (Figure  III.2)  indicate  fundamental  differences  in  the  aerosols’  microphysical 
properties  compared  with  the  results  of  other  flights  or  the  SIMMS ’93  data  which  show  less  variation  over 
the  visible  wavelengths.  The  apparent  increase  in  attenuation  at  380  nm  is  most  likely  due  to  the  presence 
of  higher  concentrations  of  small  particles  (Figure  III.3)  that  have  greater  extinction  efficiencies  as  the 
ratio  of  size-to-wavelength  (r/X)  increases  (van  de  Hulst,  1981).  Because  such  enhancements  are  most 
pronounced  in  the  results  derived  from  highly  variable  surface  measurements  we  suspect  that  tropospheric 
haze,  not  visible  against  the  "milky"  appearing  stratosphere,  may  have  contaminated  these  particular 
results.  The  sharp  decreases  evident  at  the  large-particle  end  of  the  inferred  size  spectra  are  attributed  to 
the  observed  decreases  in  opacity  for  X  >  778  nm. 

Finally,  we  speculate  that  the  apparent  longevity  of  the  relatively  large  Pinatubo  aerosols  in  the 
Arctic  may  be  due  to  the  vaporization  and  regrowth  processes  discussed  by  Hofmann  and  Rosen  (1987), 
possibly  augmented  by  general  circulation  patterns  that  favor  the  heating  of  the  Arctic  stratosphere  during 
winter  ( Robock  and  Mao,  1992)  and  the  accumulation  of  aerosols  at  high  latitudes. 
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III.2  LEAD  CHARACTERISTICS  FROM  CO-LOCATED 
AVHRR.  DMSP  OLS,  LANDSAT  TM,  AND  ERS-1  SAR 

In  Section  I  we  describe  theory  and  observations  outlining  the  effects  of  spatial  resolution  on  lead 
statistics.  Here  we  apply  and  extend  those  results  to  additional  data  sets.  To  compare  whether  substantial 
differences  might  exist  in  lead  statistics  derived  from  different  image  types,  we  have  assembled  a  data  set 
consisting  of  AVHRR  direct-readout  imagery,  Defense  Meteorological  Satellite  Program  (DMSP) 
Operational  Line-scan  System  (OLS)  data,  Landsat  Thematic  Mapper  (TM)  images,  and  ERS-1  synthetic 
aperture  radar  (SAR)  data  for  15  April,  16  April,  and  18  April  1992.  The  TM  and  SAR  imagery  were 
co-registered  as  part  of  a  separate  NASA-funded  effort.  The  AVHRR  imagery  were  obtained  from  Ron 
Lindsay  (Univ.  of  Washington),  with  additional  AVHRR  and  OLS  provided  by  Florence  Fetterer  of  NRL. 
Pixel  sizes  for  the  georeferenced  data  were  1  km  (AVHRR),  0.5  km  (OLS),  and  25  m  (SAR  and  TM). 
For  comparison,  the  AVHRR  and  OLS  sections  corresponding  the  the  SAR  and  TM  coverages  were 
subsectioned  out,  and  the  SAR  and  TM  resampled  to  1  km  pixels.  For  a  variety  of  reasons,  not  all  data 
sets  were  co-registered  with  each  other.  The  matching  data  sets  are  OLS  and  AVHRR  thermal  channels 
for  15  April;  AVHRR  (all  channels)  and  TM  (all  channels  for  16  April;  SAR  (15  April)  coregistered  with 
the  16  April  AVHRR  and  TM;  and  all  AVHRR  and  TM  channels  for  18  April.  The  image  sets  for  15 
and  16  April  are  discussed  here. 

The  lead-statistics  program  (hereafter  "LEADSTATS")  described  in  Section  1.3.1  was  used  to 
calculate  leads  data  from  each  image  section.  The  images  were  first  converted  to  binary  "lead/not-lead" 
images  using  an  interactive  thresholding  procedure.  Thresholds  were  set  to  capture  as  many  apparent  leads 
as  possible  (open  water  and  thin-ice  leads)  in  each  image.  While  this  procedure  is  subjective,  it  allowed 
us  to  maximize  the  information  content  in  each  image.  The  resulting  binary  images  were  then  processed. 

III.2.1  AVHRR  AND  OLS:  EFFECTS  OF  SPATIAL  AND  RADIOMETRIC  RESOLUTION 

Figures  III.4a  and  III.4b  show  co-located  OLS  and  AVHRR  thermal-band  imagery,  with  a  portion 
of  each  image  resampled  (nearest-neighbor  sampling)  to  equivalent  pixel  sizes  of  250  m  (Figure  III.5a  and 
III.5b).  In  this  example,  the  greater  radiometric  resolution  of  AVHRR  tends  to  outweigh  the  effects  of 
lower  spatial  resolution;  the  AVHRR  data  capture  the  same  lead  patterns  as  visible  in  the  OLS,  as  well 
as  showing  additional  lower-contrast  leads  (either  reffozen  or  narrow  open-water  leads)  that  are  not 
resolved  in  the  OLS  image.  However,  the  effect  of  the  lower  spatial  resolution  of  AVHRR  vs.  OLS  is 
shown  in  the  apparent  lead  widths  in  the  two  images,  where  the  AVHRR  data  suggests  wider  leads  than 
is  indicated  in  the  OLS.  Both  the  radiometric  and  spatial  resolution  effects  are  depicted  in  a  transect 
through  the  two  images  (Figure  III.6).  In  this  example,  most  of  the  leads  are  2  digital  numbers  (DN)  or 
less  of  the  ice  background  DN  value  in  the  OLS  data.  The  tendency  for  AVHRR  to  overestimate  lead 
width  relative  to  OLS  is  apparent  at  locations  150  and  725  along  the  transect.  This  difference  obviously 
depends  on  the  threshold  selected  to  define  leads.  In  Figure  III.7,  for  example,  selecting  a  threshold  for 
each  image  type  to  capture  only  the  warmest  pixels  yields  quite  similar  results.  Along  this  transect,  the 
thresholds  used  in  Figure  III.7  results  in  8.0%  of  the  transect  mapped  as  leads  in  the  AVHRR  data,  and 
7  0%  in  the  OLS  data  When  the  same  thresholds  are  applied  to  the  full  image  subsets  in  Figure  ili.4, 
the  difference  in  lead  area  is  greater  (total  lead  area  of  18.2%  in  the  OLS  and  28.5%  in  the  AVHRR), 
indicating  the  difficulty  of  selecting  a  single  threshold  to  yield  uniform  results  across  even  a  relatively 
small  portion  of  images. 

The  detail  in  the  AVHRR  data  due  to  the  greater  radiometric  resolution  compared  to  OLS  has  a 
large  effect  on  the  ability  of  the  LEADSTATS  routine  to  define  leads  with  high  confidence.  For  example, 
when  cloud  cover  is  excluded  from  the  images  shown  in  Figure  III.4,  and  thresholds  manually  selected 
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to  highlight  leads,  LEADSTATS  identifies  only  a  single  clearly-defined  lead  and  3  ill-determined  leads. 
In  the  AVHRR  image.  87  leads  are  clearly  defined,  with  another  130  ill-determined  leads.  Wnile  this 
difference  might  be  reduced  by  trying  a  variety  of  thresholds,  it  points  out  that  mapping  of  lead  area  using 
spatial  patterns  and  context  is  even  more  sensitive  to  radiometric  information  than  is  lead-area  mapping 
using  thresholds  only. 

III.2.2  LANDSAT  THEMATIC  MAPPER,  AVHRR,  AND  SAR:  EFFECTS  OF  RESAMPLING, 
WAVELENGTH,  AND  SCALE 

As  noted  above,  a  set  of  Landsat  TM,  AVHRR,  and  SAR  data  were  assembled  to  allow  lead 
features  to  be  intercompared  at  different  spectral  wavelengths  and  spatial  sampling  (Figure  III. 8). 
Obvious  features  of  this  set  are  the  close  correspondence  between  the  AVHRR  and  Landsat  TM  Ch.  2 
image  (originally  30  m  field-of-view  resampled  to  1  km  pixel  size)  and  the  utility  of  the  Landsat  thermal 
channel  (Ch.  6;  originally  200  m  field-of-view)  to  discriminate  open  water  from  refrozen  leads.  Also 
apparent  are  large  differences  between  the  SAR  image  and  the  TM  and  AVHRR  data.  The  SAR  image 
was  acquired  one  day  earlier,  but  the  lack  of  direct  correspondence  to  the  TM  and  AVHRR  imagery  is 
typical  of  data  compared  elsewhere.  Figure  III  .9  highlights  the  differences  between  the  SAR  and  TM  Ch. 
2.  These  differences  are  apparently  due  to  a  large  range  in  surface  backscatter  that  is  not  necessarily 
related  to  ice  thickness,  whereas  the  albedo  and  temperature  information  in  the  TM  and  AVHRR  data 
appear  to  be  more  directly  related  to  open-water  and  refrozen  leads.  The  intermediate  backscatter  of  areas 
indicated  as  open  water  in  the  TM  and  AVHRR  suggests  that  the  open-water  areas  are  being  roughened 
by  wind,  and  thus  not  spectrally  distinguishable  from  thin  ice.  Differences  between  ERS-1  SAR  and 
optical  and  thermal-wavelength  data  are  discussed  in  further  detail  in  the  following  section  on  the 
evolution  of  lead  patterns  in  response  to  atmospheric  circulation.  A  comparison  of  Landsat  thermal 
channel  6  and  AVHRR  thermal  channel  4  (Figure  III.  10)  points  out  the  effects  of  the  higher  radiometric 
resolution  of  AVHRR,  which  shows  more  lead  detail  in  spite  of  AVHRR’s  lower  spatial  resolution.  The 
effects  of  this  lower  spatial  resolution  on  apparent  lead  width  is  also  demonstrated  in  this  figure. 

While  the  TM  and  AVHRR  data  are  qualitatively  similar  when  mapped  to  equivalent  pixel  sizes, 
comparisons  of  statistics  of  lead-covered  area  can  vary  considerably  depending  on  the  threshold  used  to 
define  lead  area  in  the  AVHRR  data.  Selecting  a  threshold  that  delineates  most  of  the  leads  in  the 
AVHRR  data  yields  substantial  overestimates  of  lead  area  compared  to  lead  area  mapped  from  TM  data 
gridded  to  a  50  m  pixel  size.  For  example,  selecting  thresholds  in  this  manner  for  AVHRR  Ch.  2  and  Ch. 
4  yields  lead  areas  of  4.7%  and  3.0%,  compared  to  0.8%  in  the  50  m  TM  Ch.  2  data.  To  explore  this 
issue  of  how  lead  statistics  intercompare  across  different  fields-of-view,  we  apply  a  method  discussed  in 
Key  et  al.  (1993b)  and  in  Section  1.3. 5. 2  of  this  report.  In  this  approach,  it  was  noted  that  lead  area 
mapped  using  a  single  threshold  tends  to  decrease  roughly  linearly  with  pixel  size.  This  decrease  is  due 
to  the  progressive  loss  of  radiometric  contrast  of  lead  pixels  with  surrounding  ice  pixels  as  spatial 
resolution  decreases.  Thus,  by  calculating  lead  area  L  at  resolution  K,  and  K',  where  K’  is  a  resolution 
degraded  following  some  spatial  interpolation  method,  then  the  resulting  slope  can  be  used  to  estimate  L 
at  some  finer  spatial  resolution. 

To  test  this  approach  as  a  means  of  intercomparing  lead  area  estimated  at  difference  pixel  sizes, 
we  degraded  the  1  km  AVHRR  for  16  April  to  pixel  sizes  of  2,  3,  and  4  km  using  bilinear  interpolation. 
The  TM  data  were  also  degraded  to  50,  100,  200,  500,  and  1000  m  pixels.  In  keeping  with  the 
assumption  of  loss  of  lead  area  with  decreasing  resolution,  a  threshold  should  be  selected  such  that  only 
the  darkest  (in  ordeal  wavelengths)  or  warmest  lead  pixels  are  flagged  as  lead  area.  The  alternate 
approach  of  selecting  a  threshold  to  maximize  the  apparent  lead  coverage  produces  an  overestimate  of  lead 
area.  Given  that  radiometric  contrast  decreases  as  leads  no  longer  fill  a  single  pixel,  attempting  to  select 
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an  intermediate  threshold  value  then  places  entire  pixels  into  a  lead-area  class,  when  in  fact  only  a  portion 
of  that  pixel  (or  perhaps  none  of  the  pixel,  depending  on  the  spatial  interpolation  used)  actually  consists 
of  lead  area.  Thus,  when  as  above,  we  select  a  threshold  in  the  AVHRR  data  to  highlight  the  lead 
patterns,  the  result  was  about  5  times  more  lead  area  than  actually  exists,  based  on  the  Landsat  TM  data. 
While  this  effect  of  threshold  selection  may  be  intuitively  obvious,  the  tendency  nontheless  is  present  to 
select  thresholds  that  give  the  most  visually  pleasing  maps  of  lead  networks.  However,  selecting  a 
threshold  by  the  "minimum  rule"  given  above  is  not  only  more  realistic  of  the  true  situation,  but  also  lends 
itself  to  automated  processing,  where  the  image  can  be  scanned  for  minimum  reflectances  or  maximum 
temperatures  to  provide  the  appropriate  thresholds. 

Thresholds  selected  in  this  manner  produced  the  following  lead-area  proportions  in  the  full- 
resolution  and  degraded  AVHRR  and  Landsat  data.  In  the  case  of  the  Landsat  data,  a  threshold  was  first 
selected  using  the  50  m  image,  and  then  applied  to  the  degraded  data,  e.g.,  as  a  test  of  what  happens  when 
the  true  reflectence  of  an  open  lead  is  known. 


Table  III.  1.  Lead  area  as  a  proportion  of  total  ice  area,  estimated  using  a  "minimum-rule"  threshold  from 
AVHRR  imagery  degraded  from  1  km  to  4  km  pixel  sizes  using  bilinear  interpolation. 


Channel 

1  km 

2  km 

3  km 

4  km 

Ch.  2 

0.63% 

0.31% 

0% 

0% 

Ch.  4 

0.28% 

0.19% 

0.14% 

0% 

Table  III.  2.  Lead  area  as  a  proportion  of  total  ice  area,  estimated  using  a  "minimum-rule”  threshold  from 
Landsat  TM  imagery  degraded  from  50  m  to  1  km  pixel  sizes  using  bilinear  interpolation. 


Channel  50  m  100  m  200  m  500  m  1  km 

Landsat  TM  2  0.78%  0.72%  0.63%  0.44%  0.23% 


The  change  in  lead  area  with  pixel  size  is  essentially  linear  in  the  AVHRR  data.  From  the  slopes  of  these 
lines,  we  can  estimate  the  lead  area  at  a  pixel  size  of  50  m,  which  we  can  take  to  represent  the  true  areal 
coverage  of  leads.  Applying  this  to  the  AVHRR  Ch.  2  and  Ch.  4  data  yields  estimated  lead  area  at  50 
m  pixel  sizes  of  0.93%  for  Ch.  2  and  0.37%  for  Ch.  4.  The  combination  of  the  minimum-rule  threshold 
selection  and  this  adjustment  for  pixel  size  takes  the  lead  area  proportion  in  the  correct  direction,  although 
with  the  small  lead  fractions  in  these  data,  the  effect  of  the  change  is  fairly  small.  Using  the  Landsat  data 
as,  for  example,  a  MODIS-type  product  where  one  wished  to  estimate  lead  area  from  500  m  data,  and  was 
able  to  use  a  subset  of  a  higher-resolution  MODIS  channel  or  coincident  Landsat  data  to  define  a  "true" 
lead  reflectance  threshold,  then  by  degrading  the  500  m  data  to  1  km,  the  resulting  slope  results  in  an 
estimated  lead  area  of  0.63%  for  a  50  m  pixel  size.  This  is  clearly  an  improvement  over  simply  using  the 
lead-area  estimate  from  the  500  m  data  alone.  If  nothing  else,  this  approach  helps  to  avoid  overestimating 
lead  area  from  low-resolution  data,  while  offering  an  adjustment  that  compensates  at  least  partly  for  the 
loss  of  contrast  with  interpolation. 
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The  detail  present  in  the  full-resolution  TM  data  (not  shown)  make  the  TM  imagery  particularly 
well-suited  to  test  the  effects  on  lead  statistics  of  spatial  resampling.  Using  LEADSTATS,  177  well- 
deftned  leads  (representing  1.7%  of  the  total  area)  were  identified  in  the  full-resolution  TM  image,  43 
leads  (also  1.7%  of  the  image  area)  at  100  m  pixel  size  resampled  using  nearest-neighbor  sampling,  and 
48  leads  (1.5%  of  the  area)  at  100  m  resampled  using  bilinear  interpolation.  The  difference  between  the 
sampling  methods  is  typical  of  the  loss  of  spatial  pattern  but  retention  of  radiometric  integrity  offered  by 
nearest-neighbor  sampling,  whereas  bilinear  interpolation  tends  to  retain  the  linearity  of  leads  but  due  to 
interpolation,  a  portion  of  lead  pixels  are  lost  due  to  radiometric  averaging  with  surrounding  ice  pixels. 
This  trade-off  between  retaining  spatial  patterns  and  preserving  the  spectral  signatures  will  depend  on  lead 
width,  since  the  detectability  of  narrow  leads  will  be  more  affected  by  spatial  interpolation  than  will  large 
leads.  The  effects  of  interpolation  are  also  dependent  on  the  contrast  between  leads  and  the  surrounding 
ice  cover. 

To  investigate  the  effects  of  AVHRR  spectral  channel  on  lead  retrieval,  lead  statistics  were 
calculated  from  AVHRR  channels  1-4  at  1  km  pixel  size,  and  from  the  SAR  image  resampled  to  1  km 
pixels  (e.g.,  the  image  sets  in  Figure  III.  8).  More  lead  segments  are  defined  using  AVHRR  Ch.  1  (39 
leads)  than  Ch.  2  (29  leads)  perhaps  due  to  the  slightly  greater  contrast  between  open  water/thin  ice  and 
snowcover  in  Ch.  2.  The  total  lead  area  is  about  the  same  (approximately  4%  of  the  image  area).  With 
the  AVHRR  thermal  channel  4,  15  leads  are  found  (about  1%  of  the  area).  The  difference  between  the 
lead-area  estimates  at  the  optical  and  thermal  wavelengths  appears  to  be  a  reasonable  estimate  of  the 
differences  between  the  proportions  of  open  water  and  thin  ice  within  the  pack.  In  the  SAR  image,  the 
total  number  of  lead  fragments  identified  (15)  and  total  lead  area  (1%)  were  similar  to  the  AVHRR  data, 
although  the  actual  areas  classified  as  lead  bore  little  resemblance  to  leads  mapped  in  the  AVHRR  data. 

m.2.3  EVOLUTION  OF  LEAD  PATTERNS  IN  AVHRR  AND  ERS-1  SAR  DATA 

Spacebome  SAR  is  the  only  sensor  that  can  produce  sequential  coverage  of  the  ice  pack  under 
all  weather  conditions  at  spatial  resolutions  suitable  to  detect  individual  leads.  In  principal  then,  SAR  data 
should  be  well-suited  for  observing  the  evolution  of  leads  under  different  wind  and  ice-motion  regimes. 
However,  as  demonstrated  above,  substantial  uncertainties  exist  in  how  leads  actually  appear  in  terms  of 
backscatter  properties,  compared  to  lead  patterns  apparent  in  optical  and  thermal-wavelength  imagery. 
To  investigate  whether  the  analysis  of  a  sequence  of  SAR  images  rather  than  a  single  scene  can  aid  in 
interpreting  lead  patterns,  a  time  series  of  SAR  data  has  been  interpreted  for  locations  in  the  Beaufort  Sea 
during  periods  of  strong  synoptic  activity  during  15-31  October  1991,  and  for  an  annual  cycle  from 
October  1991  through  October  1992  (Maslanik  et  al.,  1993;  Heinrichs  et  al„  1993).  Results  from  the 
October  case  study  is  discussed  below.  This  work  was  also  supported  by  NASA  funds,  where  the 
SAR-derived  lead  fractions  and  ice  concentrations  were  compared  to  SSM/I,  AVHRR,  and  ice  model 
output.  Here,  we  focus  on  the  leads-mapping  issue. 

For  the  15-31  October  case  study,  three  sets  of  SAR  imagery  resampled  to  100  m  pixels  were 
acquired  for  three  locations  in  the  central  Beaufort  Sea.  Three  to  four  overlapping  images  were  selected 
for  each  site.  As  in  the  examples  cited  above,  a  backscatter  threshold  was  selected  for  each  image  to 
discriminate  as  well  possible  open  water/thin  ice  from  first-year  and  old  ice  in  the  uncalibrated  SAR 
images.  Statistics  on  numbers  of  leads,  lead  width,  and  orientation  were  derived  for  each  image  using  the 
LEADSTATS  program.  Visual  interpretation  of  the  SAR  time  series  sets  provides  the  best  indication  of 
how  lead  patterns  evolve  under  the  different  wind  regimes  in  the  October  case.  Ice  motion  is  slow 
enough  to  allow  observation  of  the  same  features  in  the  ice  pack  over  most  of  the  sampling  period.  At 
all  three  sites,  no  leads  are  visible  in  imagery  acquired  on  17  October,  prior  to  the  passages  of  deep 
low-pressure  systems  on  21-23  October.  As  the  lows  move  through  the  region,  leads  develop  in  the 
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central  Beaufort  oriented  in  nearly  the  same  direction  as  the  geostrophic  winds.  By  26  October,  the  lead 
orientations  have  changed  to  primarily  north-south  in  the  northern  and  central  Beaufort  sites,  with  no 
leads  visible  at  the  southern  site.  The  same  pattern  is  present  on  29  October  at  each  site. 

Using  the  LEADSTATS  routine,  the  maximum  amount  of  any  of  the  SAR  images  classified 
definitely  as  lead  area  did  not  exceed  0.7%.  When  all  pixels  with  backscatter  values  below  the 
classification  threshold  were  included,  the  estimate  of  open  water/thin  ice  area  is  still  small  for  all  the  SAR 
images  studied  for  the  October  case,  with  little  evidence  of  a  large-scale  opening  of  the  pack  due  to 
passage  of  the  strong  synoptic  systems.  Either  little  change  in  open-water  fraction  occurs  due  to  the 
variable  winds,  or  the  changing  ice  conditions  are  not  apparent  in  the  SAR  data.  The  time  span  between 
images  is  long  enough  to  allow  refreezing.  However,  since  some  leads  remain  visible  over  5  to  9  days, 
it  is  likely  that  some  evidence  of  a  substantial  opening  of  leads  would  be  apparent.  It  is  also  possible 
that  openings  occur  that  are  sufficiently  small  to  be  undetected  ai  the  SAR  pixel  size  (100  m  in  the  Lo-Res 
imagery  used  here).  Analyses  using  full-resolution  (25  m)  ERS- 1  SAR  imagery  for  a  different  application 
(Steffen  et  al.,  1993)  did  not  suggest  that  the  resampling  to  100  m  resulted  in  a  significant  change  in 
open-water/thin  ice  fraction.  However,  this  comparison  needs  to  be  repeated  for  the  specific  cases 
considered  here. 

Discrimination  of  leads  in  SAR  data  is  complicated  by  the  complex  backscatter  properties  of  open 
water  and  thin  ice.  As  noted  above,  the  main  sources  of  uncertainty  and  lack  of  unique  backscatter 
signature  are  the  effects  of  wind  speed  (affecting  open  water  areas),  and  variability  in  surface  dielectric 
properties  and  surface  roughness  that  can  occur  in  thin  ice  forming  in  leads.  Open-water  areas  can  have 
very  low  backscatter  under  calm  conditions  with  backscatter  increasing  as  wind  speed  increases,  such  that 
open-water  leads  can  be  indistinguishable  from  surrounding  ice,  as  is  apparently  the  case  in  the 
comparison  of  Landsat  TM,  AVHRR,  and  SAR  discussed  earlier.  If  wind  speeds  are  high  enough,  and 
perhaps  depending  on  wind  direction  and  lead  orientation,  open-water  areas  can  be  the  highest-backscatter 
features  in  SAR  images  (Steffen,  et  al.,  1993). 

AVHRR  data  for  the  October  case  study  provide  a  qualitative  look  at  sea  ice  conditions,  although 
cloudiness  is  too  extensive  to  permit  digital  classification  of  ice  concentration  from  AVHRR  during  much 
of  the  period.  In  general,  as  was  the  case  with  the  comparison  with  TM  and  SAR,  the  AVHRR  data 
suggest  greater  amounts  of  leads  than  the  SAR  data  indicate.  Given  the  relationships  of  SAR  backscatter 
to  a  number  of  factors  such  as  wind  speed,  surface  brine  expulsion,  and  frost-flower  formation,  it  is  quite 
possible  that  leads  apparent  in  AVHRR  thermal-channel  data  will  not  be  visible  in  the  SAR.  For  example, 
Figure  III.  1 1  shows  co-registered  AVHRR  and  SAR  coverages  for  29  October.  While  the  large  leads  in 
the  SAR  are  also  visible  in  the  AVHRR  thermal-band  data,  other  leads  apparent  in  the  AVHRR  image 
do  not  appear  in  the  SAR.  The  physical  temperatures  of  the  leads  in  this  AVHRR  image  average  about 
5  K  warmer  than  the  surrounding  ice  (245  K  vs.  240  K),  indicating  that  the  leads  either  consist  of  young 
ice  (rather  than  open  water),  or  are  narrower  open-water  leads  that,  through  resampling  and  interpolation, 
appear  as  wider,  lower-temperature  features.  A  threshold  classification  of  this  portion  of  the  AVHRR 
image  corresponding  the  the  SAR  scene  yields  about  7%  of  the  area  as  open-water/thin  ice  ("warm"  ice' 
65%  first-year  ice  (mid-  temperature  ice),  and  25%  old  ice  ("cool"  ice). 

Since  the  apparent  lead  temperature  at  the  AVHRR  field-of-view  is  a  function  of  lead  width,  lead- 
surface  temperature,  and  ice-surface  temperature  (e.g.,  the  observed  pixel  temperature  is  a  mixture  of 
temperatures  of  the  lead  and  ice  surface),  then  knowledge  of  lead  width  reduces  this  system  to  a  single 
unknown  (lead-surface  temperature)  if  we  assume  that  accurate  lead  width  can  be  estimated  from  SAR, 
and  ice-surface  temperature  can  be  reasonably  estimated  from  AVHRR  data  elsewhere  in  the  image.  Thus, 

Tl  =  (Tb  -  T,Pj)/Pl  , 
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where  TL  is  the  lead-surface  tempeature,  TB  is  the  observed  AVHRR  pixel  temperature,  T,  is  the  thick- ice 
surface  temperature.  Pj  is  the  proportion  of  thick  ice  in  the  pixel,  and  PL  is  the  proportion  of  lead  in  the 
pixel.  For  the  image  in  Figure  III.  1 1,  using  the  large  lead  in  the  center  of  the  images,  T„  =  approximately 
243  K,  T,  =  240  K,  and  since  lead  width  from  the  S AR  is  about  400  m  and  the  AVHRR  pixels  are  2  km 
by  2  km,  PL  =  400  m/2000  m  and  P,  =  (2000  m  -  400  m)/2000  m.  From  the  above  equation,  TL  thus 
equals  255  K,  suggesting  that  the  lead  is  in  fact  refrozen.  The  single  example  given  here  could  be  refined 
by  using  full-resolution  (12.5  m)  SAR  to  get  estimate  lead  width  as  accurately  as  possible,  since  if  the  lead 
is  in  fact  200  m  wide,  then  Tt  is  near  270  K.  Alternatively,  if  it  is  known  that  the  lead  is  actually  open 
water,  such  as  if  wind  roughening  yielded  a  very  high  backscatter  clearly  due  to  wavelets  on  the  water 
surface,  then  this  knowledge  can  be  used  to  constrain  the  estimate  of  lead  width.  In  the  case  given  above, 
for  example,  the  lead  would  be  unlikely  to  be  much  smaller  than  200  m,  since  T,  then  becomes  greater 
than  the  freezing  temperature.  A  100  m  lead  would  require  a  TL  of  about  300  K.  Given  this  surface 
temperature,  and  also  wind  speed,  radiation,  conductivity,  and  ocean  heat  flux,  a  rough  estimate  could  be 
made  of  the  ice  thickness  within  the  lead.  Another  application  of  this  pixel  "unmixing"  through 
combination  of  SAR  and  AVHRR  might  be  to  define  a  true  reflectance  or  surface  temperature  for  leads, 
which  could  then  be  used  in  the  procedure  discussed  earlier  to  correct  lead-area  estmates  for  pixel  size. 
This  would  be  equivalent  to  the  degraded-Landsat  example,  where  the  best  threshold  to  use  was  known 
(since  it  was  determined  from  the  50  m  TM  data).  In  this  case,  however,  the  threshold  does  not  require 
a  high-resolution  image  such  as  Landsat,  and  can  make  use  of  the  wider  availability  of  SAR  data. 

As  one  would  expect  from  the  operation  of  the  leads-mapping  routine,  area  classified  as  leads  is 
substantially  smaller  than  total  open  water/thin  ice  area  defined  by  a  simple  threshold.  In  the  case  of  the 
SAR  data  studied  here,  the  "spatial  coherence"  criteria  that  the  leads-mapping  routine  uses  to  define  leads 
is  particularly  important,  since  manual  interpretation  of  the  SAR  data  shows  that  much  of  the  area 
classified  as  open  water/thin  retain  a  consistently  low  backscatter  over  the  full  15-day  study  period,  which 
is  unlikely  given  the  expected  refreezing  of  open  water  areas.  The  lead  areas,  however,  evolve  in  a  logical 
manner  over  time.  These  results  indicate  that,  as  suggested  by  other  investigators,  some  of  the  uncertainty 
in  SAR  classification  can  be  alleviated  using  pattern  recognition  tools  in  addition  to  simple  backscatter 
thresholds. 

m.2.4  SUMMARY 

Comparisons  of  the  different  image  types  allows  us  to  assess  the  effects  of  spatial  resolution, 
radiometric  resolution,  resampling  method,  and  temporal  information  on  leads  mapping.  AVHRR  and 
Landsat  Thematic  Mapper  imagery  (resampled  to  match  the  AVHRR  pixel  size)  contain  essentially  the 
same  information  in  terms  of  discimination  of  leads.  The  choice  of  resampling  methods  (bilinear 
interpolation  or  nearest-neighbor  sampling)  depends  on  whether  retention  of  spatial  patterns  or  radiometric 
information  is  most  critical. 

While  SAR  data  clearly  contain  leads  information,  comparison  to  the  other  image  types  shows  that 
the  variable  backscatter  of  thin  ice  and  open  water  in  the  single-channel,  single-polarization  ERS-1  SAR 
substantially  complicate  the  interpretation  of  leads,  and  can  potentially  cause  large  errors  when  automated 
leads  mapping  routines.  Accurate  discrimination  of  lead  features  in  single-  channel  SAR  data  will  likely 
require  spatial  operators.  However,  even  when  spatial  features  such  as  linearity  and  orientation  are  used, 
SAR  classification  is  limited  by  a  variety  of  factors  that  can  render  lead  backscatter  indistinguishable  from 
the  backscatter  of  the  surrounding  ice.  It  is  likely  that  the  best  classifiers  will  consider  both  spatial  and 
temporal  patterns  and  changes  in  backscatter,  as  well  as  the  actual  backscatter  values  themselves. 

For  cases  where  the  same  leads  are  apparent  in  both  SAR  data  and  AVHRR,  then  it  should  be 
possible  to  estimate  the  lead  surface  temperature  from  the  AVHRR  thermal  channels,  given  the  "true"  lead 
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width  as  estimated  from  the  SAR  data.  An  example  is  given  for  a  matched  set  of  AVHRR  and  SAR 
images.  This  combination  could  unambiguously  define  open-water  leads  from  refrozen  leads.  Also,  with 
additional  information  such  as  air  temperature,  winds,  and  oceanic  heat  flux,  this  surface  temperature  could 
be  used  to  derive  the  ice  thickness  within  refrozen  leads.  In  practical  terms,  the  accuracy  of  this  approach 
is  limited  by  differences  in  acquisition  time  between  the  SAR  and  AVHRR  data  and,  to  a  much  less 
extent,  errors  in  retrievals  of  surface  temperature  from  AVHRR.  Other  potentially  valuable  approaches 
exist  that  take  advantage  of  combinations  of  data  and  ancillary  information  to,  for  example,  estimate  lead 
fraction  when  wind  speeds  are  known  (Heinrichs  et  al.,  1993),  or  to  calculate  wind  speeds  when  leads  are 
known  to  be  open  water  (Steffen  et  al.,  1993). 
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Fig.  III.5.  Enlarged  portions  of  OLS  (a,  top)  and  AVHRR  (b,  bottom)  for  15  April  1992,  with  both 
image  types  resampled  to  250  m  pixels.  Location  of  the  transect  plotted  in  Figure  III.6  is  indicated  on 
(b). 
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AVHRR,  OLS  transect  (15  Apr.) 


Transect  (250  m) 


Fig.  III.6.  Transect  of  AVHRR  and  OLS  digital  numbers  along  the  transect  shown  in  Fig.  III.5.  The 
AVHRR  transect  is  plotted  above  the  OLS  transect,  "in  offset  of  50  was  added  to  the  OLS  digital 
numbers  for  plotting  purposes. 
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AVHRR,  OLS,  lead  thresh.  (15  Apr.) 


Transect  (250  m) 


Fig.  III.7.  Comparison  of  portions  of  the  transect  mapped  as  leads  using  AVHRR  and  OLS  digital- 
number  thresholds.  Digital  numbers  greater  than  the  selected  thresholds  were  mapped  to  a  value  of  200. 
Apparent  are  the  limited  radiometric  range  of  the  OLS  data,  effect  of  threshold  on  estimated  lead  width, 
and  the  good  spatial  registration  between  the  two  images. 
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Fig.  III.  11.  Co-located  SAR  and  AVHRR  coverage  for  the  central  Beaufort  Sea.  Typical  lead  widths  in 
the  SAR  data  are  about  500  m.  Some  lead  patterns  probably  consisting  of  thin,  warmer  ice  are  apparent 
in  the  AVHRR  image  but  not  in  the  SAR  data. 
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Here  we  report  on  the  application  and  validation  of  the  1ST  retrieval  procedure  using  a  suite  of 
in  situ  data  collected  during  the  May-June  1992  Seasonal  Sea  Ice  Monitoring  and  Modelling  Site  (SIMMS) 
field  campaign.  Satellite-derived  estimates  of  1ST  are  compared  to  near-surface  air  temperatures,  surface 
temperatures  derived  from  upwelling  longwave  (broadband)  radiation  measurements,  and  from  a 
thermocouple  just  below  the  snow  surface.  Since  the  satellite  field-of-view  (FOV)  is  at  best  1  1  km, 
comparisons  are  also  made  with  surface  temperature  estimates  collected  over  an  approximately  1  km2  area 
with  a  hand-held  infrared  thermometer.  The  problems  encountered  are  not  unlike  those  experienced  in 
the  validation  of  sea  surface  temperatures,  where  there  are  a  variety  of  ways  of  measuring  the  surface 
temperature  in  situ.  Note  again  that  the  "surface"  temperature  refers  to  the  skin  or  radiating  temperature 
of  the  snow  or  ice  surface. 

III.3.1  METHODOLOGY 

Surface  microclimate  data  were  collected  during  the  SIMMS’92  (Seasonal  sea-ice  Modeling  and 
Monitoring  Site)  field  experiment.  The  first-year  ice  (FYI)  site  [74°  41.66’  N,  95°  35.22’  W]  was  the 
focus  of  the  measurement  program  and  operated  from  19- April  to  26-June.  Air  temperatures  were 
measured  with  thermocouple  sensors.  All  sensors  are  accurate  to  within  approximately  0.3°  C.  The 
sensors  were  housed  within  ventilated  phychrometer  shieldings.  While  air  temperatures  were  measured 
at  five  vertical  levels,  the  air  temperature  sensor  used  for  this  study  (Tal)  ranged  in  height  between  54 
cm  and  57  cm  above  the  snow  surface  from  May  10  to  May  24,  and  between  47  cm  and  65.5  cm  after 
May  24  to  the  end  of  the  experiment.  As  with  air  temperatures,  snow  temperatures  were  measured  in 
profile  within  the  snow  cover  extending  from  the  snow/ice  interface  to  near  the  snow  surface,  using  a 
thermocouple  epoxied  in  the  tip  of  white  plastic  tubing  (20  cm  x  0.5  cm).  All  sensors  and  leads  were 
painted  white.  The  temperature  from  the  thermocouple  nearest  the  surface,  the  depth  of  which  varied 
from  0-3  cm,  is  used  here. 

Downwelling  sky  and  surface  emitted  broadband  infrared  radiation  (4-50  pm)  were  measured  with 
Eppley  Precision  Infrared  Radiometers.  One  of  these  pyrgeometers  was  mounted  on  an  extension  arm 
approximately  8.9  m  from  the  snow  surface  on  an  instrument  tower.  A  sky  facing  pyrgeometer  was 
mounted  on  a  post,  approximately  35  m  west  of  the  instrument  tower.  Instrument  height  was 
approximately  1.5  m  above  the  snow  surface.  The  flux  measurement  was  corrected  for  the  infrared 
radiation  emitted  by  the  thermopile  using  the  Stefan  Boltzman  equation  and  a  measure  of  the  internal 
instrument  temperature  as  recorded  using  the  precision  thermistor  (YSI  44031)  housed  within  the 
instrument.  Thermistor  tolerance  is  approximately  ±0.125°  C  between  -10°  and  -20°  C;  however  the 
measurement  error  is  probably  closer  to  ±0.2°  C  when  logged  to  the  21X  micrologger.  Regardless  of 
theses  specifications,  the  accuracy  of  the  instrument  is  difficult  to  quantify  because  of  possible  heating  of 
the  instrument  dome  due  to  absorbed  incident  short-wave  radiation;  consequently,  the  temperature  of  the 
dome  was  also  monitored  in  order  to  assess  the  degree  to  which  any  heating  of  the  dome  may  bias  the 
flux  measurements. 

The  surface  thermodynamic  temperature  based  on  pyrgeometer  data,  Tpyrg,  was  estimated  from  the 
upwelling  longwave  radiation  emitted  by  the  surface,  L\,urf  according  to  the  Stefan-Boltzman  law: 
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where  Tpyrg  is  in  degrees  C,  a  is  Boltzman’s  constant,  and  e  is  the  emissivity,  taken  to  be  0.99  here. 
However,  because  the  emissivity  is  less  than  unity,  the  upwelling  radiation  measured  by  the  sensor,  Lt, 
is  the  sum  of  the  radiation  emitted  by  the  surface,  LTJurf,  and  that  portion  of  the  downwelling  atmospheric 
radiation,  Li,  that  is  reflected  by  the  surface: 

L\  •  L\UHt  +  (1-e  )Ll 

so  that  this  reflected  component  must  be  removed  from  the  upwelling  radiance  before  calculating  the 
surface  temperature.  The  discussion  in  the  appendix  suggests  that  the  uncertainty  surrounding  our  Tpyrg 
estimate  is  approximately  ±0.2°  C. 

In  order  to  characterize  the  spatial  variability  of  the  temperature  field  over  an  area  comparable  to 
an  AVI1RR  pixel,  skin  temperatures  measured  with  an  infrared  (IR)  thermometer  (T,Rlhtm)  and  snow/ice 
interface  temperatures  (Ta  i)  were  measured  along  transects  within  the  1  x  1  km  sample  site  and  the 
multiyear  site.  Each  set  of  surface  observations  at  the  FYI  site  consisted  of  measurements  spaced  200  m 
apart  along  two  randomly-selected  1  km  transects,  with  the  time  between  measurements  kept  to  a 
minimum.  The  typical  time  required  to  cover  the  ten  stops  was  about  90  minutes,  with  sampling  times 
selected  to  correspond  to  AVHRR  overpasses  whenever  possible.  Tm,herm  was  measured  with  an  Everest"” 
hand-held  IR  thermometer;  a  non-contact  instrument  that  determines  a  brightness  temperature  of  the  object 
within  the  instrument’s  field-of-view  based  on  received  radiation  in  the  8-14  pin  range.  The  instrument 
is  factory-calibrated  to  yield  a  representative  accuracy  of  0.5°C  in  an  operating  environment  of  0°  -  50°C. 
The  manufacturer  estimates  that  accuracy  below  0°C  is  approximately  the  same.  The  instrument  was 
tested  periodically  by  measuring  the  temperature  of  fresh  water  at  its  freezing  point  in  a  slush  bucket, 
where  the  IR  thermometer  typically  yielded  temperatures  of  +/-  0.3°C.  During  measurements  along  the 
transects,  the  IR  thermometer  was  allowed  to  reach  equilibrium  temperature  with  its  surroundings. 
Measurements  were  taken  with  the  IR  thermometer  held  about  1  m  from  the  surface  at  an  angle  to  the 
surface  of  approximately  45°.  Emissivity  of  the  snow-covered  sea  ice  was  set  at  0.99. 

For  the  retrieval  of  1ST  we  use  the  methodology  of  Key  and  Haefliger  (1992),  as  described  in  an 
earlier  section.  Local  Area  Coverage  (1.1  km  field-of-view  at  nadir)  data  from  the  NOAA-1 1  and  NOAA- 
12  satellites  collected  by  Atmosphere  Environment  Service  are  used.  First-order  calibration  was  performed 
following  the  methods  described  in  NOAA  (1991a).  Additional  corrections  were  applied  to  the  data  to 
account  for  the  nonlinear  response  of  the  thermal  channels  (Weinreb  et  al.,  1990;  NOAA,  1991b).  The 
AVHRR  scan  angle  ranges  from  0°  to  approximately  55°.  Since  the  1ST  retrieval  can  only  be  done  for 
clear  sky  conditions,  and  because  automated  cloud  detection  in  the  polar  regions  is  difficult  at  best,  images 
that  are  clear  over  the  first-year  ice  site  are  selected  through  a  visual  analysis  of  various  combinations  of 
the  AVHRR  visible,  near-infrared,  and  thermal  channels.  However,  it  appears  that  even  this  manual 
interpretation  of  the  imagery  may  not  be  adequate  for  detecting  low-level  ice  crystal  precipitation 
("diamond  dust")  and  very  thin  stratus. 

III.3.2  RESULTS 

The  near-surface  air  temperature,  Taln  the  surface  temperature  from  the  thermocouple  buried  just 
beneath  the  snow  surface,  Ttmv,  and  the  temperature  based  on  upwelling  longwave  radiation,  Tpyrg,  are 
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shown  in  Figure  III.  12  for  May  and  June,  at  a  solar  time  of  approximately  1000.  Values  for  T)now 
represent  all  measurents  where  the  thermocouple  was  less  than  3  cm  below  the  snow-air  interface.  The 
snow  temperature  is  typically  higher  than  the  other  two  temperature  measurements  in  the  early  part  of  the 
experiment  due  to  the  insulating  effect  of  the  overlying  snow.  Snow  temperatures  are  sometimes  greater 
than  the  melting  point  of  the  snow  after  the  initial  stages  of  melt  onset  in  late  June  due  to  warming  of  the 
sensor  by  increased  transmission  of  solar  radiation  throughout  the  snow  layer.  Table  III.3  provides  a 
statistical  comparison  of  the  different  measures  of  "surface"  temperature.  The  two  periods  in  May  reflect 
early  spring  (8-18  May)  and  the  transition  to  late  spring  conditions  (20-24  May),  when  surface  air 
temperatures  become  warmer  than  T, ,,  signifying  a  change  in  the  direction  of  heat  transport  toward 
warming  of  the  ice  surface  at  the  snow/ice  interface.  TIRlhtnn  and  Tt  i  were  grouped,  and  averaged  over  half 
hour  intervals  to  coincide  with  the  half  hour  averages  of  T  and  Tair 


Table  III.3.  Comparison  of  mean  temperatures  (°C)  and  range  of  differences  for  8-25  May  at  the  FYI  site. 


Time  Period 

8-18  May 

20-25  May 

8-25  May 

Mean  Temperatures: 

-13.39 

Tair  -14.74 

-6.74 

-7.72 

-11.76 

-13.02 

-13.90 

-6.52 

-12.09 

Tsi  -8.92 

-6.40 

-8.15 

Temperature  Differences: 

T  -  T 

M  IR therm  ‘  pyrg • 

Mean  Difference  -0.51 

0.22 

-0.33 

Standard  Dev.  0.62 

0.62 

0.69 

T  -  T  * 

Mean  Difference  0.84 

1.20 

0.93 

Standard  Dev.  1.32 

0.44 

1.17 

T  -  T  : 

pyg  air 

Mean  Difference  1.35 

0.98 

1.26 

Standard  Dev.  1.19 

0.90 

1.13 

T  -  T  • 

1  IRtherm  1  si * 

Mean  Difference  -5.69 

-0.42 

-4.08 

Standard  Dev.  2.26 

2.33 

3.36 
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Values  of  TAVHRR  are  plotted  with  Tpyrg  in  Figure  III.  13  for  May  and  June.  The  temperature 
reported  as  TAVHRR  is  a  mean  of  10  pixels  around  the  FYI  site  whose  fields-of-view  contained  first-year  ice 
only  (i.e.,  no  land),  covering  approximately  13  km2,  depending  on  the  sensor  scan  angle.  Gaps  in  the  time 
series  result  from  periods  of  extensive  cloud  cover.  The  maximum  time  difference  between  the  in  situ 
and  satellite  observations  is  approximately  20  minutes.  Differences  between  the  temperature  pairs  range 
from  less  than  0.1  K  to  more  than  3  K,  with  the  satellite  estimates  almost  always  less  than  the  in  situ 
values.  The  differences  between  TAVHRR  and  Tair  (not  shown)  are  less,  but  the  sign  of  the  difference  is 
usually  the  same. 

Temperatures  obtained  from  the  IR  thermometer,  TIRlherm,  are  compared  to  Tpyrg  and  the  snow-ice 
interface  temperature,  Ts  i,  in  Figure  III.  14.  The  systematic  difference  between  TIRlherm  and  Tpyrg  over  the 
8-25  May  period  is  small  (0.33°  C  difference);  Tmkfrm  measuring  slightly  lower  in  early  spring  (-0.51°  C) 
and  slightly  higher  during  the  late  spring  conditions  (0.22°  C).  The  discrepancy  is  less  than  that  shown 
in  Figure  III.  13  between  the  T  and  TAVHRR  which  is  not  entirely  unexpected  since  TAVHRR  represents  a 
sample  in  time  and  Tpyrg  a  time  average.  The  difference  between  T,Rlherm  and  Tpyrg  at  any  one  time  may  also 
be  large,  as  shown  by  both  the  fairly  large  standard  deviations  (Table  III.3).  Absolute  differences  ranged 
between  0.02°  C  and  1.80°  C.  T3i  follows  the  general  increase  in  heat  input  to  the  surface  energy  budget, 
with  less  day  to  day  variability  because  of  the  insulating  affect  of  the  overlying  snow  cover  (Figure  III.  14). 
Snow  depth  averaged  26±9  cm  from  8-25  May. 

III.3.3  DISCUSSION 

The  temperature  differences  observed  in  our  study  are  likely  due  to  undetected  clouds  in  the 
imagery,  the  spatial  and  temporal  variability  of  the  temperature  field,  incorrect  assumptions  concerning 
the  surface  emissivity  and  atmospheric  conditions,  inaccuracies  in  the  model  used  in  the  development  of 
the  satellite  retrieval  algorithm,  and  instrument  error. 

III.3.3.1  Undetected  Clouds 

Despite  the  multi-spectral  approach  to  the  manual  selection  of  clear  images,  there  are  conditions 
where  condensed  water  simply  cannot  be  detected  with  the  spectral  information  that  the  AVHRR  provides. 
This  is  particularly  true  for  very  thin  clouds  such  as  Arctic  stratus.  Even  more  difficult  to  detect  are  low- 
level  water  or  ice  fogs,  the  most  common  being  ice  crystal  precipitation  in  winter  and  early  spring.  The 
problem  in  their  detection  is  that  they  usually  exist  within  the  low-level  temperature  inversion  and  may 
result  in  top-of-atmosphere  radiances  very  close  to  what  would  be  observed  in  their  absence;  i.e.,  their 
radiative  properties,  both  shortwave  and  longwave,  are  similar  to  those  of  the  surface.  This  often  equates 
to  a  top-of-atmosphere  temperature  difference  of  a  few  degrees  or  less. 

A  re-examination  of  the  satellite  data  after  the  initial  manual  cloud  cover  assessment  and  a 
comparison  with  cloud  observations  taken  at  the  meteorological  station  Resolute  Airport  indicate  that  some 
form  of  thin  cloud  may  actually  have  been  present  over  the  FYI  site  at  the  times  of  the  AVHRR 
acquisitions.  However,  as  cloud  conditions  may  differ  over  only  a  few  kilometers,  it  is  impossible  to  state 
conclusively  how  often  this  problem  influences  our  analysis.  While  cloud  type  and  opacity  were  estimated 
for  the  entire  sky  hemisphere  during  the  IR  thermometer  measurements,  additional  information  is  needed 
concerning  clouds  that  lie  along  the  path  between  the  satellite  and  the  surface  temperature  measurements. 
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111.3.3.2  Spatial  and  Diurnal  Variability  of  Surface  Temperature 

Variability  of  Tmhtrm  within  the  1  km  site  arises  primarily  from  differences  in  snow  depth,  the 
shadowing  effects  of  sastrugi-like  features,  and  changes  in  temperature  during  the  time  required  to  cover 
the  transects.  Shown  are  ail  observations  of  TIRlllfrm  throughout  the  period,  plotted  by  local  time.  The 
apparent  relationship  between  solar  zenith  angle  and  skin  temperature  is  biased  somewhat  by  sampling  date 
for  times  between  1600  hrs  and  2000  hrs  since  these  observations  were  all  collected  on  a  single  day.  The 
greatest  source  of  spatial  variability  is  likely  due  to  shadowing.  Measurements  taken  with  the  shaded 
versus  sunlit  aspects  of  sastrugi  (typical  height  of  10  -  15  cm)  filling  the  entire  field  of  view  of  the  IR 
thermometer  yielded  temperature  differences  of  as  much  as  10°C  in  clear-sky  conditions,  but  only  1°C 
when  overcast.  In  the  normal  sampling,  the  IR  thermometer  was  held  further  away  from  the  surface  to 
yield  a  larger  field-of-view  that  typically  encompassed  a  mixture  of  shadow  and  sunlit  areas,  but  not 
necessarily  in  equal  proportions. 

Figure  III.  15  shows  the  spatial  variability  of  TAVHRR  in  June  over  a  10  pixel  area  around  the  FYI 
site.  There  were  no  open  water  or  thin  ice  (significantly  less  than  1  m)  areas  in  the  scene.  The  standard 
deviation  in  the  temperatures  measured  with  the  IR  thermometer  is  up  to  twice  that  of  the  AVHRR -derived 
ISTs.  Minimizing  the  time  to  sample  the  transects  would  reduce  the  variability  of  TlR,Ktrm. 

The  surface  temperature  may  change  rapidly  in  response  to  radiative  forcing.  T  represents  an 
average  of  this  temperature  over  a  half  hour  period.  While  TIRllunn  also  represents  an  average,  its  value 
may  be  biased  by  the  non- systematic  number  of  samples  within  the  averaging  period  (ranging  from  1  to 
6).  the  irregular  time  increments  between  samples,  and  by  spatial  variability. 

111.3.3.3  Assumptions  Made  in  the  Retrieval  Procedures 

A  number  of  questions  arise  concerning  the  assumed  properties  of  the  surface  and  atmosphere 
affecting  the  AVHRR  data  and,  potentially,  the  pyrgeometer  measurements.  One  such  assumption 
concerns  the  emissivity  of  the  surface  for  the  IR  thermometer  measurements.  All  reported  observations 
of  Tuuk'rm  are  based  on  an  emissivity  setting  of  0.990,  which  is  the  highest  setting  possible  on  the 
instrument  This  value  is  probably  appropriate  for  the  broadband  emissivity  of  snow  with  grain  sizes  in 
the  100  -  200  pm  range,  but  not  for  the  8-14  pm  spectral  range  of  the  IR  thermometer.  In  that  range,  the 
emissivity  of  snow  is  closer  to  0.993  ( Dozier  and  Warren,  1982).  The  relationship  between 
thermodynamic  tempertaure  T  (K),  emissivity  e,  and  radiance  L  (milliwatts/nr-steradian-cm  ')  is  derived 
from  the  Planck  function  as 


where  c,  and  c2  are  constants  (c,  =  1. 1910659  x  10 5  milliwatts/m^steradian-cm^  and  c,  =  1.438833  cm-K) 
and  X  is  the  wave  number  (cm1).  While  the  difference  between  the  assumed  and  probable  snow 
emissivity  is  not  large,  it  does  result  in  a  temperature  overestimate  of  approximately  0.3°C. 

Additionally,  the  tipwelling  radiance  measured  by  the  instrument  is  the  sum  of  the  radiation 
emitted  by  the  surface  and  the  downwelling  atmospheric  radiation  reflected  by  the  surface.  (We  consider 
any  depletion  or  addition  of  radiation  from  the  atmosphere  itself  to  be  insignificant  given  the  short 
atmospheric  path  from  the  sensor  to  the  surface.)  While  the  amount  of  downwelling  atmospheric  radiation 
reflected  by  the  surface  is  relatively  small  given  the  high  surface  emissivity,  it  does  increase  the  radiance 
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received  by  the  sensor.  If  corrections  for  this  effect  were  made,  values  of  TIRlhtrm  would  be  even  lower 
and  closer  to  the  satellite-derived  temperatures  by  an  estimated  0.1  K  to  0.2  K. 

Assumptions  concerning  atmospheric  conditions  could  play  a  role  in  explaining  the  differences 
between  the  satellite-derived  and  in  situ  temperatures.  As  described  earlier,  the  retrieval  algorithm  is 
based  on  ice  station  data  from  the  central  Arctic  Ocean.  The  coefficients  were  derived  for  three  seasons 
and  in  some  sense  represent  the  mean  conditions  during  those  periods.  How  do  these  mean  conditions 
compare  with  those  observed  during  SIMMS ’92?  May  falls  into  the  transition  season  of  Key  and 
Haefliger  (1992),  while  June  falls  into  their  summer  season.  Differences  between  the  ice  station  and  the 
SIMMS  location  (represented  by  nearby  Resolute)  are  shown  in  Table  III .4  in  for  total  precipitable  water, 
near-surface  air  temperature  and  aerosol  optical  depths.  No  information  is  available  on  the  actual  aerosol 
optical  depth  at  the  drifting  ice  station  so  that  the  value  shown  in  the  table  is  the  assumed  amount  used 
in  model  calculations.  While  the  surface  air  temperatures  at  the  two  locations  (and  times)  are  different, 
the  range  of  temperatures  used  in  the  development  of  the  1ST  retrieval  algorithm  is  similar  to  that 
experienced  during  SIMMS ’92.  The  water  vapor  amounts  are  also  similar.  Aerosol  amounts  are  different, 
primarily  as  a  result  of  stratospheric  aerosols  from  the  eruption  of  Mt.  Pinatubo,  as  observed  during 
AG  ASP  IV  (Arctic  Gas  and  Aerosol  Sampling  Program)  over  the  Beaufort  Sea,  1992  ( Stone  et  al.,  1993 
and  Section  III.  1  of  this  report).  Although  the  overall  effect  of  aerosols  on  the  attenuation  of  upwelling 
longwave  radiation  is  small,  an  aerosol  amount  greater  than  that  expected  would  produce  a  lower  top-of- 
atmosphere  radiance  which,  especially  at  significantly  off-nadir  views,  would  result  in  an  underestimate 
of  1ST. 


Table  III.4.  Mean  precipitable  water,  surface  temperatures,  and  aerosol  optical  depths  at  a  drifting  ice 
station  and  Resolute,  N.W.T. 


Ice  Station 
Transition1 

Summer2 

Resolute 

May 

June 

Precipitable  Water 

4.6 

6.6 

4.2 

7.1 

(mm) 

Surface  Temperature 

-19.4 

-0.9 

-13.8 

-2.7 

(C) 

Aerosol  Optical  Depth 

rx \  — r~~rz — r~r. 

0.073 

: — i  - 

0.07 

2.5 

2.0 

'April,  May,  September,  1987 

2June  -  August,  1987 

3Assumed  value  for  model  calculations 


Lastly,  there  may  be  a  bias  in  the  retrieved  ISTs  due  to  AVHRR  calibration  errors.  This  may  not 
be  trivial;  for  example  the  nonlinear  calibration  alone  can  make  a  difference  of  more  than  1°C.  Overall, 
however,  we  expect  calibration  errors  to  be  much  less,  on  the  order  of  0.1  K  to  0.2  K. 

III.3.3.4  Possible  Errors  in  the  Calculation  of  Tpyrg 

Given  that  the  adjustments  to  TIRllurn  and  TAVHRR  just  discussed  will  decrease  the  differences 
between  the  two  measurements,  and  given  that  both  of  these  temperatures  are  typically  lower  than  that 
derived  from  the  average  upwelling  longwave  broadband  radiation  Tpyrg,  we  are  left  wondering  if  there  are 
any  reasons  why  Tpyrg  would  be  biased. 
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The  vertical  placement  of  the  pyrgeometer  at  8.9  m  above  the  surface  could  play  a  role.  Early 
in  the  season  when  temperature  inversions  are  common,  the  atmosphere  between  the  surface  and  the  sensor 
is  warmer  than  that  of  the  surface.  Given  that  the  atmosphere  i-  relatively  warm,  that  the  surface 
emissivity  is  less  than  unity,  and  that  relative  humidities  near  the  surface  are  high,  one  would  expect  an 
increase  in  longwave  radiation  received  at  the  sensor  over  that  emitted  by  the  surface.  Radiative  transfer 
calculations  of  upwelling  longwave  flux  confirm  this,  although  the  increase  in  the  estimated  physical 
temperature  of  the  surface  is  small,  0.3°  -  0.4°  C  for  a  surface  temperature  of  -13°C.  Additionally,  the 
large  field  of  view  associated  with  the  pyrgeometer  samples  the  emitted  infrared  radiation  from  a  much 
larger  area  than  the  hand-held  pyrometer,  and  the  ER  emitted  by  the  support  tower  infra-structure  and  other 
’non-snow’  objects  are  also  sampled. 

Our  derivation  of  L  using  the  pyrgeometer  assumes  that  any  solar  heating  of  the  instrument  dome, 
and  subsequent  re-radiation  to  the  thermopile  is  negligible,  and  that  the  thermopile  temperature  is  closely 
approximated  by  the  case  temperature.  There  is  little  reason  to  believe  a  large  discrepancy  in  temperature 
between  the  thermopile  and  case  thermistor  given  their  close  proximity.  The  net  effect  of  the  solar  heating 
of  the  dome  is  to  increase  the  calculated  value  of  L.  Solar  heating  of  the  dome  is  not  expected  for  the 
inverted  pyrgeometer  (LT),  and  in  fact  the  case  and  dome  temperature  differed  on  average  by  only  0.10° 
C  between  8-25  May,  with  the  dome  being  actually  cooler  than  the  case.  The  opposite  was  observed  for 
the  sky  facing  pyrgeometer.  Slight  heating  of  the  dome  for  Li  was  observed,  but  is  not  considered  to  be 
a  factor  that  significantly  influences  T  rg. 

III.3.4  CONCLUSIONS 

In  an  effort  to  validate  the  satellite  retrievals  of  ice  surface  temperature,  differences  between 
AVHRR-derived  ice  surface  temperature.  TAVHKK,  and  the  radiating  temperature  derived  from  measurements 
of  upwelling  longwave  (broadband)  radiation,  Tpyrg,  during  May  and  June  in  the  Canadian  Arctic  were 
observed  to  range  from  less  than  0. 1°C  to  more  than  3°C  with  TAVHKK  always  less  than  Tpyrg.  Similarly,  the 
mean  temperatures  of  spatially-distributed  measurements  made  by  a  hand-held  IR  thermometer,  TIRlhfrm, 
were  typically  less  than  Tpyrg,  although  the  differences  were  not  as  great  as  between  TAVHRR  and  Tpyrg.  The 
temperature  measured  by  a  thermocouple  placed  approximately  1  cm  below  the  snow-air  interface 
illustrates  the  insulating  quality  of  the  overlying  snow,  being  higher  than  the  radiating  temperature  in  the 
early  part  of  the  season. 

While  these  results  can  be  explained  to  a  limited  extent  by  instrument  calibration,  incorrect 
assumptions  of  the  surface  emissivity  and  atmospheric  conditions,  and  model  inaccuracies,  the  main  issue 
with  a  validation  exercise  such  as  this  is  in  the  definition  of  the  "correct"  surface  temperature  and  of  the 
method  chosen  to  measure  this  temperature.  The  temperature  of  interest  is  the  radiating  temperature,  not 
the  near-surface  snow  temperature  or  the  near-surface  air  temperature.  In  theory  the  upwelling  longwave 
radiation  can  be  used  to  determine  this  temperature,  but  in  practice  care  must  be  taken  to  obtain  an 
accurate  value.  The  spatially-averaged  temperatures  measured  with  the  IR  thermometer  suggest  a  possible 
negative  bias  in  the  AVHRR-derived  temperatures,  with  a  possible  simple  correction  for  this  bias. 
However,  the  spatial  and  temporal  variations  in  skin  temperatures  observed  during  the  sampling  of  the 
1  km  transects  with  the  IR  thermometer  are  too  great  to  permit  a  determination  of  a  bias  with  any 
certainty.  We  conclude  that  there  is  probably  a  negative  bias  in  TAVHKK  as  computed  here,  on  the  order  of 
0.5°  C  to  1°  C.  To  determine  the  actual  magnitude  and  source  of  this  bias  we  need  a  well-calibrated, 
airborne  radiometer  measuring  in  the  same  spectral  bands  as  the  AVHRR  (to  reduce  the  atmophseric 
effects  and  to  obtain  adequate  spatial  coverage),  combined  with  refined  surface  observations  that  include 
faster  sampling,  more  precise  instrumentation,  and  more  detailed  observations  on  cloud  properties  and 
distributions  relative  to  the  AVHRR  scan. 


Part  III:  Analysis  of  LEADEX  and  Other  Data 


118 


Surface  Skin,  Air,  and  Snow  Temps 


Fig.  HI.  12.  Near-surface  air  temperature,  temperature  measured  by  a  thermocouple  just 
below  the  surface,  and  temperature  derived  from  upwelling  longwave  radiation  at  the  FYI 
site. 
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AVHRR  and  In  Situ  Surface  Temperatures 


Date 


Fig.  HI.  13.  Surface  temperatures  estimated  from  upwelling  longwave  radiation 
measured  at  the  surface  and  from  the  AVHRR  thermal  channels  at  the  Fust-year 
ice  site  during  May  and  June. 
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In  Situ  Surface  Temperatures 


Date 


Fig.  HI.  14.  Surface  temperatures  measured  with  the  IR  thermometer  (spatial  means) 
and  those  based  on  upwelling  longwave  radiation.  Also  shown  are  snow-ice  interface 
temperatures  measured  (means). 
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Standard  Deviation  (C) 
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AVHRR  Skin  Temperature  Variation 


Fig.  HI.  15.  Variability  of  satellite-derived  surface  temperature  within  a  10-pixel  area 
around  the  FYI  site. 
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The  method  described  in  an  earlier  section  for  the  retrieval  of  lead  widths  and  orientations  has 
one  serious  shortcoming:  it  is  based  on  a  linear  regression  procedure  and  is  therefore  inapplicable  to 
leads  or  lead  fragments  that  are  curved.  Additionally,  the  origin  of  the  coordinate  system  upon  which 
correlations  and  regression  fits  were  based  was  arbitrary,  but  linear  correlation  is  not  insensitive  to 
rotation.  We  have  therefore  investigated  the  use  of  a  recently-developed  method  for  the  retrieval  of 
lead  geometrical  characteristics:  skeletal  modeling  ( Banfield ,  1992). 

Briefly,  this  procedure  uses  skeletons,  a  tool  from  mathematical  morphology,  to  provide  a 
useful  structural  description  of  the  shape  of  a  lead.  The  skeleton  of  an  object  is  the  set  of  points  that 
describes  the  "center"  of  the  object,  in  this  case  the  centers  of  the  largest  open  balls  that  touch  the 
boundary  of  the  lead  at  two  or  more  locations. 

The  procedure  was  applied  to  the  AVHRR  image  shown  in  Figure  HI.  16.  Only  clear  sky  ice 
areas  were  examined.  Cloud  clearing  was  done  by  applying  a  threshold  to  the  channel  3  (3.7  pm) 
reflectances,  which  was  estimated  by  subtracting  the  Pianck  radiance  for  that  channel  based  on  the 
temperature  measured  in  channel  4  (11  pm)  and  dividing  the  remaining  radiance  by  the  solar  radiance 
in  that  band. 

A  binary  image  (lead/ice)  was  then  used  as  input  to  the  code  provided  by  J.  Banfield.  The 
results  are  shown  in  Figure  III.  17  for  lead  area,  lead  width,  and  lead  skeletal  length  (the  length  along 
the  lead  as  opposed  to  the  length  of  the  main  diagonal). 
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Fig.  m.17.  Lead  areas,  widths,  and  skeletal  lengths  for  the  image  in  Figure  III.  16,  as  determined  by 
the  skeletal  modeling  method. 
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Part  iv:  Summary  of  Accomplishments 


IV.  1  ACCOMPLISHMENTS 

Over  the  course  of  this  three-year  project,  empirical  studies  concerning  scale  relationships  in  the 
retrieval  of  sea  ice  lead  statistics  have  been  performed,  as  have  modeling  investigations  of  atmospheric 
influences  on  the  satellite  signal.  Additionally,  statistical  models  that  describe  the  scaling  properties  of 
leads  have  been  developed.  The  empirical  studies  have  been  based  on  Landsat  imagery,  while  the 
atmospheric  models  have  been  specific  to  the  AVHRR.  Submarine  sonar  data  have  been  used  in  the 
statistical  model  development.  Specific  accomplishments  include: 

(1)  The  parameterizations  of  clouds,  haze,  and  atmospheric  chemical  constituents  in  the  LOWTRAN-7 
radiative  transfer  model  have  been  reviewed.  Atmospheric  temperature  and  humidity  profiles  for 
the  Arctic  have  been  constructed  from  Soviet  ice  island  data  and  were  used  in  the  construction 
of  three-season  "standard"  atmospheres  for  the  central  Arctic. 

(2)  Resampling  methods  have  been  tested  on  simulated,  AVHRR,  and  Landsat  images,  as  have  the 
effects  of  digital  enhancements.  Nearest-neighbor  resampling  has  been  shown  to  be  the  most 
effective  in  maintaining  the  geometrical  and  spectral  characteristics  of  leads. 

(3)  Empirical  relationships  between  pixel  size  and  lead  width  have  been  illustrated. 

(4)  Procedures  for  the  retrieval  of  lead  statistics  have  been  developed  and  applied  to  Landsat  imagery 
successively  degraded  to  more  coarse  resolutions.  Applications  to  ERS-1  SAR  data  have  also 
been  explored,  although  a  single-channel  SAR  does  not  appear  to  provide  adequate  spectral 
information  for  ice  type  discrimination.  New  methods  of  retrieving  lead  statistics  utilizing 
mathematical  morphological  techniques  have  been  applied  to  AVHRR  data  collected  during 
LEADEX. 

(5)  The  relationship  between  "apparent"  lead  widths  measured  along  a  transect  (e.g.,  from  submarine 
sonar  or  as  a  sampling  method  for  satellite  imagery)  and  the  "true”  lead  width  distribution  has 
been  formalized  in  a  statistical  sense,  so  that  one  distribution  may  be  obtained  from  the  other. 
Submarine  sonar  data  have  been  analyzed  in  this  context. 

(6)  A  statistical  model  has  been  developed  for  the  retrieval  of  lead  area  fraction  from  measurements 
along  a  line;  e.g.,  a  submarine  sonar  transect  or  a  lineal  sampling  method  for  satellite  images. 

(7)  The  effects  of  atmosphere/surface  conditions  on  the  AVHRR- measured  radiance  in  the  thermal 
channels  have  been  examined  in  terms  of  the  thermal  contrast  between  leads  and  the  surrounding 
ice  pack.  Surfaces  include  open  water,  5  cm,  15  cm,  and  2  m  thick  ice.  Atmospheric  conditions 
include  clear  sky  with  haze,  cirrus,  and  low-level  ice  crystal  plumes.  The  relationship  between 
atmospheric  optical  depth  and  lead  size  has  been  quantified. 

(8)  Information  on  Arctic  aerosol  optical  depths  during  LEADEX  has  been  collected  by  two  of  the 
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investigators  on  the  NOAA  P-3  with  sun  photometers.  This  data  has  been  used  to  assess  the 
extent  of  tropospheric  and  stratospheric  aerosols  and  their  effect  on  lead  detection  from  satellite. 

(9)  Satellite  and  in  situ  data  collected  during  LEADEX  and  SIMMS ’92  have  been  used  to  test  the 
theoretical  and  empirical  models/methods  developed  during  the  first  two  project  years,  and  shows 
these  models/methods  to  be  generally  valid. 

Additionally,  two  workshops  for  the  satellite  remote  sensing  investigators  of  the  Leads  ARI  prior 
to  LEADEX  were  hosted  by  this  group  in  Boulder.  Seven  referreed  papers  have  been  published,  with 
three  others  submitted  for  publication.  Two  graduate  students  have  been  supported  part-time  over  the 
course  of  the  project. 


IV.2  RECOMMENDED  RESEARCH 

Future  work  should  include  an  investigation  of  the  shortwave  channels  of  the  AVHRR  and  OLS 
sensors,  following  a  methodology  similar  to  that  presented  here  for  the  thermal  channels;  i.e.,  the  effects 
of  aerosols,  clouds  and  ice  crystal  precipitation  on  upwelling  radiation  should  be  examined.  Concerning 
the  retrieval  of  lead  geometrical  characteristics,  using  the  visible  and  near-infrared  channels  should  aid  in 
the  unmixing  of  pixels  containing  a  lead  of  unknown  ice  thickness  and  "background"  ice,  where  having 
measurements  at  a  variety  of  wavelengths  should  reduce  the  number  of  possible  ice  thickness/lead  width 
combinations,  as  suggested  earlier. 
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